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Abstract 

Shell models of hydrodynamic turbulence originated in the seventies. Their main aim was to describe the statistics 
of homogeneous and isotropic turbulence in spectral space, using a simple set of ordinary differential equations. In 
. the eighties, shell models of magnetohydrodynamic (MHD) turbulence emerged based on the same principles as their 
hydrodynamic counter-part but also incorporating interactions between magnetic and velocity fields. In recent years, 
significant improvements have been made such as the inclusion of non-local interactions and appropriate definitions 
for helicities. Though shell models cannot account for the spatial complexity of MHD turbulence, their dynamics 
are not over simplified and do reflect those of real MHD turbulence including intermittency or chaotic reversals of 
large-scale modes. Furthermore, these models use realistic values for dimensionless parameters (high kinetic and 
'■ magnetic Reynolds numbers, low or high magnetic Prandtl number) allowing extended inertial range and accurate 
dissipation rate. Using modern computers it is difficult to attain an inertial range of three decades with direct numerical 
simulations, whereas eight are possible using shell models. 

In this review we set up a general mathematical framework allowing the description of any MHD shell model. 
The variety of the latter, with their advantages and weaknesses, is introduced. Finally we consider a number of 
applications, dealing with free-decaying MHD turbulence, dynamo action, Alfven waves and the Hall effect. 
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Variables u, b, p, z ± , a, u ± , b* 

External magnetic field, forcing, vorticity bo, f, oi 

Wave numbers k, p, q 

Energies, helicities, enstrophy, squared magnetic potential E", E b , H u , H b , H c , S, A 

Complex conjugation, complex conjugate of z c.c, z* 

Dimensionless numbers Re, Rm, Pm, Rh 

Viscosity, diffusivity, density v, rj, p 

Injection rate of energy, kinetic helicity, cross helicity, magnetic helicity e, £ , x, £ 

Viscous, Joule dissipation rate e v , e n 

Rotation rate, Alfven wave velocity Q., bo 

Characteristic time scales t^i, ?n, tA, t v , t, ; 

Characteristic length scales If, l v , l, ; 

Characteristic velocity and magnetic fluctuations at scale / U\, b[ 

Velocity and magnetic structure functions, scaling exponent S p (l),S b (l), £ p 

Specific wave number modulus kf, k v , k n , k ± , k\\ 

Shell common ratio A 

Kinematic growthrate T^n 

Shell variables U„, B„, Z*, £/*, B% 

Shell wave number, forcing k n , F n 

Quadratic functions Q(X), W(X, Y), W(X, Y) 

Quadratic quantities in shell n E%, E%, H„, H%, H„ A„ 

Energy transfer Tf Y , Tf Y 

Energy flux nf < , nf < , nf > 

Mode-to-mode energy transfer S (i\j\k) 

Non-locality parameter y 



Table 1: Table of notations. 



1. Introduction 



In astrophysical objects m ost fluids are electrically conductin g and generally exhibit highly turbulent motion due to 
the large dimensions involved (ISchekochihin and Cowlev[|2007l) . Such magnetohydro dynamic (MHD) turbulence is at 
the heart of the dynamo a ction generating magnetic fields in planets, stars and galaxies ( Brandenburg and Subramaniani 



2005 



2001 



Tobias et"alll201 lb . Dynamo action has been the object of several experiments dGailitis et al 1 l2000l: I S tie glitz and Miillei 

2007l:lspence etall 120071: iNataf et "alll2008tlFrick et alll2010b and~ 
is suspected to occur in nuclear reactors cooled with liquid so dium ( Plunian et all Il999t) . MHD turbulence is also 



responsible for the propagation of Alfven waves dAlfveriLll942h in presence of an external m agnetic field as in e.g. the 
solar wind . Such waves can be reproduced in MHD experiments (lAlboussiere et all 1201 II) and measured in plasma 
tokamaks ( IGekelman , Il999h . Complementary to observation and experi ment, direct numeric al simulations aimed at 



reproducing the finest details of MHD turbulence have been performed ( Miiller et all 120031) . However, one serious 



difficulty faced by simulations is that the processes i nvolved are s trongly non-linear implying, for example, that en- 
ergy is transferred over an extended range of scales ( VermaL 2004 ). This range is several orders of magnitude larger 
than what is attainable with present day or, indeed, projected computers. In this respect shell models are of primary 
importance in building up our understanding. First introduced to deal with hydrodynamic (HD) turbulence, shell 
models have now been generalized to MHD, leading to interlocked progresses of both types of model. We will now 
su mmarize the evolu tion of these ideas. 

IObukhovl(ll97ll) introduced a multilevel system of non-linear triplets to mimic the energy transport, in the spirit 
of the Richardson-Kol mogorov scenario for the energy cascade in HP turbulence. This idea has been success fully 



developed by his tea m ( GledzerJ, 11973c iDesnianskii and NovikovL 1 19741 : iGlukhovskiil 11975c iGledzer et all 119811) . At 



the same time lLorenzl ( 1 1 97 lb started from the full Fourier representation of the Navier-Stokes equations. Aiming at 
studying the statistical properties of turbulent flow with limited computer facilities, he reduced the set of equations 
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to a "very low order model". Though both approaches were different. iLorena (119721) and iGledzerl ( 1 1 9731) eventually 
derived the same shell model of HD turbulenceQ. They both applied the conservation of kinetic energy and enstrophy 
with the description of atmospheric turbulence statistics in mind. We note that these two quantities, kinetic energy and 
enstrophy, are positive definite and so are rather straightforward to define in the framework o f shell models, explain ing 



why shell models of 2D-turbulence were preferred at that time. It took a further 20 years ( Kad anoff et al.L 1 1 9951) to 



identify and adequately describe kinetic helicity (not sign definite) in shell models, leading to 3D-turbulence modeling 
(for which kinetic helicity, instead of enstrophy, is ideally conserved). 

Other low order models of turbulence appearing in the seventies were all based on the same principle: the division 
of isotropic spectral space into a set of concentric shells using only one variable per shell to characterize velocity 
fluctuations. The main difference between the m odels were the degree of locality between two interacting shells, each 



shell interacting either with one first neighbor (Ob ukhovL 1197 It iDesnianskii and NovikovL Il974t iBell and Nelkin , 



1978 ; Kerr and Siggial 19781) - with only one quadratic invariant (ki netic energy) - or two first neighbors ( Lorenz , 



197alGledzeiT 1973uGlukhovskiilll975l : |Gledzer and Makarovl 1 19791) - with two quadratic nvariants (kinetic energy 



and enstrophy). 

In the next decade, Ziminl ( 1981 ) introduced the so-called "hierarchical model of turbulence" with self-similar 
functions localized in both physical and Fourier space 0. Projecting the Navier-Stokes equations on this base of 
functions, he obtained a set of ordi nary differen tial equations organized in a hierarchical tree. By reducing this 
hierarchical tree to one vertical chain, Frik ( 1 983b Fl constructed a shell model for 2D-tur bulence - with two quadratic 
invaria nts (kinetic energy and enstrophy) - including not only local interactions as in lLorenzl (119721) and iGledzei 
(119731) but also non-local interactions. 

From the very first numerical simulations of the shell model equations, it was clear that the Kolmogorov's solution 
(or Kraichnan's in 2D) gave an unstable fixed point, and that a Kolmogorov spectrum of energy could be obtained only 
by averaging over time, as expected in real turbulence. However, such models failed to show any chaotic behavior . 



The l ink to intermittency was still missing, until the first MHD shell models were derived dFrikLll984 : Gloaguen et al 



19851) . Indeed, by doubling the degrees of freedom (adding the magnetic field), chaotic behavio ur was obtained. A 



similar effect had also been observed with temperature in shell models o f convective turbulence ([Frik, Il986lll987l) . 
Applying this idea to HD turbulence, lYamada and Ohkitanil d 1987L 1 1988b used a velocity with two real components 
instead of only one, and obtained solutions showing chaotic behavior. With such complex velocity they found inter- 
mittency statistics in excellent agreement with real HD turbulence. 



In the following years such models of HD turbulence aroused wide interest dPisarenko et all 1 1993c ICarbone , 



1994bt iBiferale et all fl99l iKadanoff et al 



ties (a three - shell p eriodicity) identified by 



dL'vov et al 
1994al) 



1995b iF rick et alj, 119951) . A spurious regularity in the spectral proper- 



Biferald dl993l was corrected either by using a slightly di fferent model 



19981) or by considering a velocity with three real components per shell instead of two dAurell et al 



Afte r the identification of kinetic helicity bv IKadanoff et all (1 19951). the first model of 3D MHD turbulence was 
derived (IBrandenburg et al. . Il996c iBasu et al.Ul998tlFrick and Sokolofn. il 998) with total energy, magnetic and cross 
helicities as quadratic invariants. Mean while, new shell models for HD turbulence were elaborated in which the 



velocity is projected onto helical modes dZimin and Hussainlll995tlBenzi et al. . Il996bl) . In such helical models the 



helicity is not correlated with the kinetic energy, cont rary to the other models. This gives rise to im portant differences 
when dealing with kinetic helicity in HD turbulence ( Stepanov et al. . 2009t Lessinnes et al. . 2011 ). It is also suitable 
to study magnetic and cross helicity in MHD turbulence ( Frick and Stepanovll2010l) . 

Within the last ten years more complex shell models have been elaborated to acco unt for characteristics peculia r to 
MHD turbulence. These models inc lude non-local int eractions, dir ectly within triads ( Plunian and Stepanov , 2007b or 



with t he help of multi-scale models ( IFrick et al. 1 12006b . anisotropy (Ni gro et alll2004l) . and the Hall effect ( IFrick et al 



2003). 



Our review is organized as follows. After a short description of MHD turbulence in Sec.|2l a general framework is 
introduced in Sec.[3]providing a description of the various shell models derived so far. This is followed in Sec.|4]by a 
review of results obtained for different applications. For the sake of clarity, HD shell models will often be introduced 



'First denoted "cascade" or "scalar" models, such models have been called "shell" in the beginning of t he nineties. 

2 In terms of contemporary scientific langage, these functions would be called wavelets, as discussed by F rick and Z iminl )1993|) . 

3 In papers, translated from Russian journals, Peter Frick was spelt as Frik P.G. and we keep each time the spelling from the cited paper. 
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before MHD shell mo dels. However, we will focus on MHD applicati o ns only. For HP applications the reader can 
refer to reviews by e.g. iBohr et all dl998l) . lBiferaie1 d2003l) . lFrickl d2003l) . lDitlevsenl d201 lb . 



2. MHD turbulence 



In this section we simply review some background infor mations necessary to address th e next secti ons. For deeper 
informations the reader can refer to reference books on HD dFriscM 1 1 995t iLesieurl 1 1 9971) and MHD dMoreau . 1 19901: 
Davidsoni 1200 It iBiskampl 120031) turbulence. 



2.1. Physical space 
2.1.1. MHD equations 

The incompressible MHD equations that govern the time evolution of the velocity u and magnetic induction b are 



(<9, -vV 2 )u = -(u ■ V)u + (b ■ V)b - Vp + f, V ■ u = 0, 
(<9,-^V 2 )b = -(u • V)b + (b ■ V)u, V ■ b = 0, 



(1) 
(2) 



where v is the viscosity, rj the magnetic diffusivity, p the total pressure (including the magnetic pressure) and f the 
flow forcing, normalized by the fluid density p. These equations are derived from the Navier-Stokes equations supple- 
mented by the Lorentz force, and Maxwell equations for which advantage has been taken of fast charge redistribution 
commonly found in liquid metals. The magnetic induction has been normalized by ^4-np such that b is given in units 
of velocity. Here we are interested only in fluctuation, assuming that u and b average to zero in space and time. 

The non-linear terms in the r.h.s. (right hand side) redistribute kinetic and magnetic energies among the full range 
of scales from the largest, defined by the system boundaries, to the smallest at which the total energy dissipates. 
Different kinds of helicities are also transfered. Such transfers are called direct or inverse, depending on whether they 
occur towards smaller or larger scales. They are also described as local or non-local depending on whether they occur 
between neighbouring scales or not. 

We speak of forced or decaying turbulence depending on whether f is different from or equal to zero and dynamo 
action when the magnetic energy does not decay in time, meaning that the transfer from kinetic to magnetic energy is 
sufficient to compensate for the magnetic dissipation. 

In presence of an external magnetic field bo (having a velocity dimension as it is normalized by ^4np), Eqs. (QJ2} 
can be rewritten by replacing b by b + bo- Introducing the so-called Elsasser variables defined as 



z ± = u ± b, 

the MHD equations become 

d,x* + (z* • V)z ± = +(b • V)z ± - Vp + r + VV + r _ vV + f=, 



(3) 



(4) 



where r* = \{v + rj) and bo is assumed to be independent of space coordinates. Provided bo is sufficiently strong, 
Eq. dH) can be linearized and thus becomes a wave equation the solutions of which are the so-called Alfven waves 
dAlfvenlll942l) . This set of equations is only symmetric for P m = 1. Taking Pm = 1 has the effect of suppressing the 
reflection of Alfven waves at the walls dSchaeffer et all 1201 II) . 



2.1.2. Quadratic invariants 

In MHD three integral quantities play a special role: the total energy, the cross helicity and the magnetic helicity. 
The total energy, E, is the sum of the kinetic energy E" and the magnetic energy E b , 



E = E U 



E\ 



7" = I f u 2 dV, E h = i f b 2 dV, 
2Jv 2Jv 



(5) 
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where V is the volume of integration. The cross helicity H c and magnetic helicity H b are defined as 



H c = I u • b dV, H b = f a • b dV, (6) 



v 

where a is the vector potential satisfying V x a = b. The absolute value of cross helicity is maximal if u and b are 
aligned (and zero if they are perpendicular). 

The quadratic quantities E, H c and H b are called invariant because in the ideal limit of a non-viscous and non- 
diffusive fluid v = n = 0, and in the absence of forcing and an external magnetic field (f = bo = 0), they are conserved 
in time, 

d,E = d,H c = d,H b = 0. (7) 

The first two conservation laws can be shown from Eqs. (Q]|2]i provided appropriate boundary conditions are used 
while the third conservation law is obtained from Eq. (0. We can also show that the first two conservation laws are 
equivalent using the following property 

C y (x-V)yrfV = 0, (8) 

which is satisfied for any divergence-free vectors x and y. Then the conservation of total energy and cross helicity 
take the following forms 

fu-(b-V)brfV + fb-(b-V)u£/V = (9) 

Jv Jv 



\ b-(u-V)urfV + f u (u V)b£/V = 0. 

Jv Jv 



(10) 



Exchanging u and b does not change E and H c . However, Eqs. (0 and (fTOt are exchanged, showing the equivalence 
of the two conservation laws. 

We note that the kinetic helicity, defined by 



codV, (11) 

v 



H u = fu 

Jv 

where co — V x u is the vorticity, is not conserved in MHD. Hence 

d,H" = 2 f b • V x (b x co) dV, (12) 
Jv 

indicating that kinetic helicity can be created or suppressed by the interaction between the vorticity and the magnetic 
field. For b = the kinetic energy and helicity are conserved. In the presence of an external magnetic field bo, the 
magnetic helicity is not conserved anymore. 

Finally, in 2D HD turbulence the kinetic helicity is always zero. Instead enstrophy 



E= \L ( ° 2 



dV (13) 



is conserved along with the kinetic energy. 

In 2D MHD turbulence magnetic helicity is always constant. Instead the square potential 



Uv'' 



2 

is conserved together with the total energy and cross helicity. 



A = - I a 2 dV (14) 
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2.1.3. Dimensionless parameters 
The Reynolds number is defined as 



Re = uil/v, 



(15) 

in Eq. (Q}, Re shows how strong the 



where u\ is a characteristic velocity and / is a characteristic scale. Putting b 
non-linear interactions are, compared with viscous dissipation. 
In MHD the magnetic Reynolds number is defined as 

Rm = uiI/t]. 

From Eq. <J3J> Rm shows how strong the non-linear interactions are, compared with magnetic dissipation. 
The ratio of both numbers is called the magnetic Prandtl number 

Pm = Rm/Re = v/77 



(16) 



(17) 



and depends only on the fluid properties. The dimensionless form of Eqs. (Q]|2]i are obtained by replacing v and 77 by 
Re 1 and RrrT 1 respectively. 

Fully developed turbulence implies high Re (greater than 10 3 at the largest scale). Dynamo action requires Rm > 1 
at the scale where the magnetic energy grows. This corresponds to typical Rm of 10 to 10 3 when Rm is calculated 
at the largest scale of the system. The latter condition is difficult to meet e xperimentally with liquid metals. Indeed 
the power necessary to run a liquid metal experiment increases as Rm 3 (IPetrelis and Fauvei l200ll) . demanding a 
considerable effort to reach Rm > 10 at the largest scale. Liquid sodium is usually used for its high conductivity, and 
for its density about unity. With Pm * 10~ 5 , reaching Rm * 10 2 would require Re =s 10 7 . 

In Fig.[TJa few typical objects are placed on the map (Re,Rm). Among them liquid metal experiments, fast breeder 
reactors, the Earth's core, Jupiter's core and the Sun's convective zone correspond to Pm ~ 10 5 to 10 6 . We note that 
such low values for Pm and also realistic Re value s remain beyond the limits of current direct numerical simulations 
dSakuraba and Roberts! l2009t lUritskv et al.Ll201C)b . 



2.1.4. Homogeneity and isotropy 

Two assumptions are usually made in order to obtain theoretical predictions for both HD and MHD turbulence. 
The first assumes homogeneity, meaning that the statistical quantities derived from the flow and magnetic fields are 
invariant under translation in physical space. This assumption fails to predict, for example, the effect of boundary 
layers. The second assumes isotropy, meaning that the statistical quantities are independent of direction. In principle 
isotropy is broken as soon as a sufficiently strong external magnetic field or rotation is applied. 

Most 3D shell models are based upon this double assumption of homogeneity and isotropy. However, several 
models have been developed to relax the assumption of isotropy in the context of Alfven waves (see Sec 14. 4. 21 . 



2.1.5. Isotropic phenomenology 



In his paper lKolmogorovl(ll941l) introduced the structure function for the velocity field 



S U P {D = (\ Ul f) 

where «/ is the longitudinal velocity increment 

Ui = [u(x + 1) - u(x)] • l/l, 



(18) 



(19) 



and ( ) denotes ensemble averaging. In HD turbulence, the power e which drives the flow at forci ng scale If, is 
transferred towards sm aller scales and, in a stationary state, is equal to the energy dissipation rate e v (IKolmogorov 
194 It lObukhovt 119411) . This direct energy cascade occurs within the inertial range corresponding to scales I such 



that l v < I < If, where l v is the viscous scale below which viscous dissipation dominates. In such an inertial range, 
assuming isotropy, a simple dimensional analysis leads to the estimation 



Uj oc (e0 



1/3 



(20) 



Figure 1: Map of typical objects in the plane (Re,Rm). DNS stands for Direct Numerical Simulation and FBR for Fast Breeder Reactor (cooled 
with liquid sodium). Yellow dashed lines are Pm isolines. 



This corresponds to an energy transfer rate uf It^i K £, where tuL = l/ui is the eddy turn-overtime. Then for any p, 

S u p (l) cc (elf 13 . (21) 

The viscous scale l v is estimated by assuming that the power e, which is injected into the fluid at some forcing scale, 
is subsequently dissipated by viscosity at scale Z v , corresponding toes; vuj //;. This leads to 

l v cc e - 1/4 v 3/4 . (22) 

In MHD turbulence at high Rm, the estimate of / ;; depends on Pm. For Pm < 1 applying the same type of phenomenol- 
ogy as above for u and b, we estimate the viscous and ohmic scales to be 

iv0c f ;"V' 4 , (23) 

where e,, and e, ; are the fractions of e that correspond to viscous and ohmic dissipation respectively (e v + e n = e) . The 
ratio of these two scales is 

/ le \ 1/4 

fioc (±) Pm-3/4. (24) 
l v \ € nl 
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For Pm > 1, assuming that magnetic energy is produced by the velocity shear which is maximum at scale l Y , the scale 
at which magnetic energy dissipates is given by l~ l u(l v ) oc rjl n 2 , leading to 



oc Pm 



-1/2 



(25) 



This immediately shows the difference between low-Pm MHD turbulence, as in liquid metal, high-Pm MHD 
turbulence, as in interstellar medium, and Pm « 1, as in direct numerical calculations. For Pm <K 1 magnetic energy 
dissipation occurs at a much larger scale than that of the kinetic energy, and vice-versa for Pm » 1. For Pm « 1, an 
example is given in Figf2] both kinetic and magnetic energies dissipate at about the same scale. 

The ratio e v /e v , about which little is known, depends on the level of magnetic energy compared to the level of the 
kinetic energy. In general we have e v » e n in experiments a nd e v < e n in real astrophysi cal objects. In MHD shell 



models, e v /e n ~ 1/10 for Pm « 10~ 5 , leading to L// v * 10 35 dPlunian and StepanovM2010l) . 



The structure functions for the magnetic field are defined similarly to Eq. (TT~8T >. as 

S b p (D = (\b,f). 



(26) 



Assuming a Kolmogorov scaling law, given by Eq. ( f20b . for both magnetic and velocity fields leads to the same scaling 



for both structure functions S" p oc S p 



In the presence of an e xternal magnetic field bp, a differen t mechanism of energy transfer occurs due to the 
interaction of Alfven waves ( Iroshnikov , 1964 : Kraichnanl 1965 ). Indeed, such an applied field leads to an additional 
time scale t& = l/bo- Provided bo is sufficiently strong, t& can be shorter than the eddy turnover time l/u/. Then energy 
transfer occurs on the Alfven time scale, leading to 



ui oc b, oc (eb Q ) 1/4 l l/ \ 



(27) 



and on to 



S u Jl)^S b Jl)^(eboy l4 l" /4 . 



(28) 



Deviation from isotropy leads to other time scales and scaling laws (see Sec. 12.2.31 ). 



2.1.6. Intennittency 

If structure functions obtained from HD experimental measurements exhibit clear scaling laws, their slopes, how- 
ever, clearly deviate from p/3 as p is increased. This is interprete d as the signa t ure of intermittent events, like bursts, 
which are not captured by the self-similarity assumption of the iKolmogorovl (119411) theory. Such intermittency is 
quantified by the scaling exponent £ p such that 

S p W*fr. (29) 

Various mode ls of intermittency have been proposed aiming at an analytical formula for the scaling exponent C, p 
( FrischL 1995 ). Here we draw attention to the elegant parameter-free formula of She and Leveque ( 1994 ) 



■HD 



= p/9 + 2(1 - (2/3y /3 ), 



(30) 



which gives ^ D = 1, consistent with the Kolmogorov "4/5" law. It is based on 

(i) the Kolmogorov refined similarity hypothesis: (uf ) oc (ef )l p ^ where e/ is the energy dissipation averaged over a 
volume Z 3 , 

(ii) log-Poisson statistics for the dissipation rate fluctuations, 

(iii) one-dimensional (filament-like) form of the ultimate dissipative structures. 

The scaling exponent given by Eq. ( 1301 is in excellent agreement with experimental measurements of isotropic 
HD turbulence. There is, however, good reason to expect that Eq. ( l30l l is not valid in MHD turbulence, even if both 
magnetic and velocity fields satisfy the same Kolmogorov scaling law given by Eq. < f2Qb . The difference comes from 
the different nature of the ultimate dissipative structures, which might be sheet-like rather than filament-like. Thus, 
after replacing hypothesis (iii) of She-Leveque by 

(iv) the ultimate dissipative structure is two-dimensional (sheet-like), 
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Horburv and Balogh 



d 1997b and 



Muller and Biskampl J2000l) proposed a MHD version of Eq. ([30} for Sf(I) 



gMHD _ pig + i_ (l/3)P/ 3 ) (31) 

which again gives c^ HD = 1 . 

When applied to Alfven wave turbulence not only must (iii) be replaced by (iv), but (i) must also be replaced by 
the Iroshnikov-Kraichan relation 

(v) < M f> ex (tf ) cc ((eibnY'y P' 4 . 

Subsequentlv lGrauer et al. I (fl994l) and lPolitano and Pouquetl dl995h developed the Alfvenic version of Eq. ( 1301 



-p/8 + l-(l/2V' /4 , (32) 



with^ = 1. 

In experiments and/or numerical simulations, accurate measurements of scaling exponents are needed in order 
to discriminate between the three formula given above by Eqs. ( 1301132b . When dealing with high order structure 
functions, which is necessary for discrimination, it is even hard to identify the appropriate range of scale s which can 



be used for the determination of the scaling laws. In this respect, significant progress has been made bv lBenzi et al 



(I1993bl) who discovered the concept of Extended Self-Similarity (ESS) while calculating high order structure functions 



from wind tunnel experimental results. The ESS takes advantage of the Kolmogorov "4/5" law 

(u]) = -^d, (33) 

so providing an exact linear relation between the third order structure function (it?) and the scale I within the inertial 
range. Instead of plotting the structure function {u P t ) versus /, leading to a power scaling f p , they plotted (m?) versus 
(m 3 ). As expected they found that the power scaled as {u p ,) oc (m 3 )^^ 3 . Furthermore, they discovered that this scaling 
held for a range of s cales I much larger (both in small and large scale directions) than that for which {u p .) oc ftp holds. 



Benz i et al. (1993b) therefore claimed an extended self-similarity range of scales. In addition, they found that by 



using this method the accuracy in the estimate of the scaling exponents was much improved. ESS was tested and 
used for the measurem ent of scaling exponents in a variety of turbulent flow conditions, including MHD turbulence 



(Rowlands et al., 2005) 



2.2. Fourier space 
2.2.1. Triads 

Assuming triply periodic boundary conditions in a cube of volume L 3 , both fields, flow and induction, can be 
expanded into discrete Fourier series: 



u(x) = u(k)e ,k x , b(x) = ^ b(k)e ,kx , k € ^Z 3 (34) 



2n. 

b(k)e"™, k e 

k k 

where u(k) and b(k) are the complex Fourier coefficients defined as 

1 r 

,-/'k x j„3 



u(k) = i J u(x)e- ik x dx 3 , b(k) = i j b(x)e _ *' x dx 3 



(35) 



The conditions u(-k) = u*(k) and b(-k) = b*(k) must be satisfied in order that u(x) and b(x) be real vectors. In 
Fourier space the divergence-free form of both u(x) and b(x) are given by 

k • u(k) = 0, k b(k) = 0, (36) 

indicating that both fields are perpendicular to k. Similarly, both Navier-Stokes and induction equations can be 
projected on to a plane perpendicular to k in Fourier space, making it possible to remove the pressure terms without 



ll 



loss of generality. These equations are ( Biskampi 12003 ) 



(d, + vk 2 )u(k) = -iP(k)- J] ((u*(q)-k)u*(p)-(b*(q)-k)b*(p)) + P(k)-f(k), 



pq 

k+p+q=0 



(3, + ^ 2 )b(k) = -iP(k)- Yj ((u*(q)-k)b*(p)-(b*(q)-k)u*(p)) 



(37) 



(38) 



pq 

k+p+q=0 



where P(k) is the operator defined by the matrix Py = <5,y — U- and corresponds to the projection of x on the plane 
perpendicular to k. Only the subset of wave numbers k, p and q satisfying k + p + q = 0, interact together. Such a 
triad is illustrated in Fig. [5] 

It can be shown that all the quadratic invariants introduced above are also conserved within each triad. Hence 
we can define energy and helicity transfer only between three modes belonging to the s ame triad. Th e formalism 
for mode-to-mode energy transfer in MHD turbulence has been developed in detail by IVermal (120041) and can be 
generalized to helicity (cross or magnetic) transfer. This formalism will be transposed to shell models in Sec. [3] 



2.2.2. Spectra 

The Fourier spectra of the quadratic quantities introduced in Sec. 12.1.21 are 



E u (k) = 



E"(k) 



H b (k) 



H c (k) 
H"(k) 



iu(k) ■ u*(k), 
1 

2 

1 / i 



b(k) b*(k), 
i(^(kxb(k))-b*(k) + c.c), 



-(u(k)-b*(k) + c.c.), 



1 



■ (i(k x u(k)) • u*(k) + c.c.) . 



(39) 
(40) 
(41) 
(42) 
(43) 



where c.c. means the complex conjugate. Their power density spectra are defined in their integral form 

X(k)= f X(k)dk, 

J\k\=k 



or discrete form 



X(k„) 



^ X(k), 



*„<|k|<*„ +1 



(44) 



(45) 



where X denotes any of the above quadratic quantities. In addition the following conditions are satisfied (IFrisch et al 
19751) 

. i n 

(46) 



\H"(k)\ < kE"(k), \H b (k)\ < k- l E b (k), \H c (k)\ < (E"(k)E b (k)) 



1/2 



Of course all these quantities depend strongly on time. However, we can look for states which are statistically station- 
ary or at least with a time dependency much larger than that of the time scale of the fluctuations (e.g. in free-decaying 
turbulence). Not only energies but also helici ties are expected to cascade. From absolute equilibrium distributions, 
which depend only on the quadratic invariants. IFrisch et alj ( 11975b found that the cascade is direct for the total energy 
and cross helicity, and in verse for the magn etic helicity. This makes a striking difference with the HD case for which, 
using the same method, iKraichnanl d 19731) found that both the kinetic energy and helicity cascades are direct. The 
results also depend on whether the MHD turbulence is forced or freely decaying , with or without the presence of an 
external magnetic field bo- 

An example of forced MHD dynamo turbulence for bo = is given in Fig. [2] From dimensional arguments, if the 
velocity and magnetic field increments satisfy m; oc b/ oc l a , then the corresponding spectral energy densities satisfy 
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Figur e 2: Plots of kinetic and magnetic energy density spectra E"(k) and E b (k), for Pm=l, obtained from a 512 3 DNS, for forced MHD turbulence. 
From lCarati et alj fcOOel) . 



E"(k) oc E b (k) oc r 2 " -1 . Assuming u/ oc bt oc (el) 1 ^, this leads to the famous "-5/3" Kolmogorov scaling law 

oc E b (k) oc ^k' 5 ' 3 , (47) 

for both kinetic and magnetic energy density spe ctra. HP turbulence experiments clearly demo nstrate the existence of 
such a scaling law over more than t hree decades (Saddoug hi and Veeravallilll994tlPope , l2000l) . This is also observed 



in DNS over more than one decade dGotoh e t all 2002J). In MHD turbulence there is, however, not such a clear inertial 
range for the magnetic energy density spectrum, as depicted in Fig. [2] Even the kinetic energy inertial range is rather 
short, less than one decade in Fig. [2] making it difficult to identify a clear scaling law. Short spectra are due to limited 
numerical resolution. Presumably future higher resolution will give rise to wider inertial range. The shape of the 
magnetic spectrum al so depends on forcing. In particular if the forcing is helical the spectrum can be peaked at the 
largest possible scale ( Brandenburg! 2001 , 20091) . Such large-scale magneti c field generation by the small-scale MHD 



turbule nce corresponds to the so-ca lled g-effect (Krause and Radlerlll980l) . For free-decaying MHD turbulence and 



Pm = 1 Miiller and Biskampl ( 2000l) and Miiller and Grappinl d2005l) found clear Kolmogorov scaling laws 



2.2.3. Spectra in presence of an external magnetic field 

As mentioned in Sec. 12.1.51 the presence of an external magnetic field bo changes the energy transfer which 
occurs at the Alfven time scale £4 rather than at the edd y turn-overtime scale Inl, provided b p = |bo| is strong enough. 
Assuming isotropy, the following spectra are expected ( Iroshnikov , 1964 : Kraichnarl 1965 ) 



E u (k) oc E b (k) oc (& e) 1/2 r 3/2 . 



(48) 



Unfortunately, in MHD turbulence experiments ( lOdier et al.L 1 1998c lAlemanv et al. , l2000l: iBavliss et all l2007h the 
values of Rm which are possible to achieve are too low (less than 10 at the largest scale) to observe a sufficiently 
wide magn etic inertial range. In t he left panel of Fig. [3] a typical magnetic energy density spectrum is plotted versus 
frequency dBourgoin et al. . 2002 ). Two slopes are observed, -1 and -11/3, neither of which can be attributed to an 
inertial range. The first slope is not easy to understand. On the other hand the second slope can be justified as follows 
dGolitsvnL fl960l: iMofFatti 1 1 96 II) . First we have to assume that the Taylor hypothesis applies in order to interpret the 
frequency as a wave number. The induction equation ©, replacing b by b + bo, where bo is taken to be independent 
of any spatial coordinate, and assuming |b| <K bp (low Rm), implies 



(d,-/7V 2 )b = (bo-V)u. 



(49) 
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Figure 3: Magnetic ener gy spectra of MHD tur bulence embedded in an external magnetic field. Left panel: spectrum obtained from the VKS 
experiment (ada pted from | Bo ur^oin et al] l2002t) 'l. The dashed line corresponds to the motor frequency (8 Hz). Right panel: Composite solar wind 
spectrum. From lBruno et al.ll200^) ! 



In a stationary statistic state the induction term (bo ■ V)u is balanced by the dissipation term rj^b, implying boku(k) 
T]k 2 b(k). Now assuming that the turbulent velocity obeys the Kolmogorov scaling law ( 1201 ). we find 



(50) 



Applying Taylor hypothesis, Eq. ( l50l l implies that E b (k) oc /~ n / 3 . We note that such line of argument assumes 
non-local interactions between the applied field bo and the small-scale turbulence. 

If the IK scaling law (l48l ) cannot be tested against experiments, at least it can be compared with observations. The 
right panel of Fig. [3] shows the magnetic energy density spectrum measured in the solar wind. The corresponding 
MHD turbulence subjected to the strong magnetic field bo emanating from the Sun, shows three successive slopes, 
again -1 at large scales, -5/3 in the inertial range, and -2.9 at the smallest scales. An interesting point here is that 
an inertial Kolmogo rov scaling law f~ 5 ^ 3 is o btained instead of th e IK scaling law /~ 3 ^ 2 (for a discussion of the two 



other slopes see e.g. IVerdini et al.l d2012al) and lHowes et alj (1201 II) ). 



In fact anisotropy plays a crucial role in Alfven wave turbulence dGoldreich an d Sridhar], 119951) . leading to mod- 
ified definitions for both the Alfven time t& oc {k\\b^) 1 and the eddy turn-over time t^i K {k±b(k ± ))~ l , where the 
subscript || and ± denote the directions parallel and perpendicular to the applied field bo- Two regimes are possible, 
depending whether t& <K fjvx or Ia ~ fjvL. They are denote d weak and strong tur bulent regime respectively. In the weak 
regime, on the basis of resonant three-wave interactions, iGaltier et al. (2000) found a cascade restric ted to the per- 
pendicular plane, with E h (k ± ) oc k] 2 . This has been numerically confirmed jBoldvrev and Perez , l2009h . In the strong 
regime the cascade occurs in both perpendicular and parallel directions. Provided the critical balance ?a = frvz. is satis- 
fied, t he magnetic energy spectrum is now expected to satisfy E h (k ± ) oc k^ 5 ^ and E h (k \i ) oc k , 2 ( Goldreich and Sridhai , 
1995 ). This seems to be well supported by solar wind measurements ( Horburv et all 20081). but still lacks numerical 



confir mation. Instead simulations give an energy spectrum E b (k ± ) oc kT^ 2 dMiiller and Grappinl 120051 : iMason et al 
20081) . It has been sug gested that such discrepancy is due to the dominance of the one Elssasser variable on the 
other (Boldvrev 



2006). However, recent results based on shell models in the perpendicular direction (Sec. I4.4l2l 



manage to reproduce the t ransition between weak and strong turbulence for a ratio varying from to 1 

( Verdini and Grappinl 2012 ). 
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2.2.4. Transfer functions 

From Eqs. ( I37H38I ) and following lVermal d2004l) . the time evolution of the Fourier modes of the kinetic and mag- 
netic energies ^"(k) and E h (k) are given by 



{d, + 2vk 2 )E"{k) = \ Yj (5""(k|p,q) + 5" 6 (k|p,q)) + !R{f(k)-u*(k)} (51) 

k+p+q=0 

(d l + 2i ] k 2 )E»(k) = I V (s ta (k|p,q)+S M (k|p,q)), (52) 
£ p,q 

k+p+q=0 

where each 5 AV (k|p, q) term represents the energy transfer rate from the modes p and q of field y, into the mode k of 
field x. They are defined as 

S MM (k|p,q) = 5""(k|p|q) + 5""(k|q|p) (53) 

r"(k|p,q) = 5""(k|p|q) + 5" i (k|q|p) (54) 

5 fc "(k|p,q) = S ta (k|p|q)+S te (k|q|p) (55) 

5^(k|p,q) = S^Cklplq) + 5**(k|q|p), (56) 

with 

S""(k|p|q) = -3{[k-u(q)][u(k)-u(p)]} (57) 

S"*Ck|plq) = +3{[k-b(q)][u(k)-b(p)]} (58) 

5 fa (k|p|q) = +3{[k-b(q)][b(k)-u(p)]} (59) 

S**(k|p|q) = -3{[k-u(q)][b(k)-b(p)]}, (60) 

and where the terms 5 vy (k|q|p) are obtained from 5 vy (k|p|q) by exchanging p and q. Each 5 " (k|p|q) term represents 
the mode-to-mode energy transfer rate from the mode p of field y into the mode k of field x, with the mode q acting 
as a mediator. 

An other way to write Eqs. d5TH52b is 

(d, + 2vk 2 )E"(k) = r""(k) + T ub (k) + 91 {f(k) ■ u*(k)} (61) 
(d, + 2?]k 2 )E h (k) = T bu (k) + T bb (k), (62) 

where the quantities T " (k) are interpreted as the energy transfer rate from all modes of the y-field into the k mode of 
the x-field. They are defined as 

7"*(k) = i V 5 fV (k|p,q). (63) 
L pq 

k+p+q=0 

The Fourier space is divided into spherical shells a„ that contain all wave vectors k such that k„ < \k\ < k n+ \. The 
energy transfer rate from the shell m of field y, to the shell n of field x is given by 

tZ= X 5Xy ( k IPl«i)- ( 64 > 

ke«„,pe«„ ( 
k+p+q=0 

Defining the kinetic and magnetic energies in shell n as 

El = I J) |u(k)| 2 , < = \Y, |b(k)|2 ' (65) 

keo„ keo„ 
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we obtain 



d,K = X( r - + 7 '-)- D " + F « 

HI 



(66) 
(67) 



with 

D|; = vJV|u(k)| 2 , D^ij^k^t, /^ = £&{f(k)-u*(k)}. (68) 

k£fl„ keo„ kea„ 

Similar definitions of shell-to-shell energy transfers were used in Mininni et al ] d2005h and Alexakis et al 1 d2005l) . A 
schematic representation of the shell-to-shell energy transfers is given in Fig. [4] (left), together with in Fig . [4] (right) 
the results obtained from the same DNS used to produce the energy spectrum of Fig. [2] The two top figures in Fig. [4] 
(right) show that the u-io-u and b-Xo-b transfers are local and forward. On the other hand the two bottom figures show 
that the b-to-u and u-to-b transfers are non-local, with a strong contribution coming from the velocity forcing shell to 
all magnetic shells. We stress that in DNS the sequence of k„ is chosen to be arithmetic in contrast to shell models in 
which the sequence is geometric. 
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Figure 4: Left panel: Schematic representation of the shell-to-s hell energy transfer s between velocity and magnetic fields. Right panel: Shell-to- 
shell energy transfers T*j„ from y in shell m to x in shell n. From lCarati et all )2006l) . 

Extending the previous definitions to fluxes, we can separate Fourier space into two parts, inside and outside a 
sphere of radius k„. We define four fluxes from y to x, from the inside / outside of the y-sphere to the inside / outside 
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of the x-sphere, 



n;T v< - 2 S- tv (k|p[q) (69) 



|k|<A-„,|p|<A„ 
k+p+q=0 



|k|>i„.|p|<*„ 
k+p+q=0 



TC y> = Yj 5 " V (k|p[q) (71) 



\k\<k,„\ P \>k„ 

k+p+q=0 



n;r v> = 2 5" v (k|plq)- (72) 



|k|>*„,|pl>*„ 
k+p+q=0 



We note that spherical she ll transfers are rather unsuited to anisotropic turbulence, e.g. with a strong bo ( ITeaca et al 



2009). lAlexakis et al.l (12007ft introduced cylindrical shells concentric with the direction of bo, and plane layers per- 



pendicular t o the latter, l eading to transfer maps between cylindrical shells on the one hand and parallel planes on the 
other hand. ITeaca et al.1 (120091) introduced ring-to-ring transfers by dividing each spherical shell into rings. This has 
the advantage of showing how the transfers are distributed with the angle between the ring and the direction of bo. 

Helicity transfer can also be defined in a way similar to that used for energy transfer. Starting from Eqs. ( 13711381 . 
the time evolution of the cross helicity H c {k) and magnetic helicity H*(k) are 



(5, + (v + ^ 2 )// c (k) = 1 J] S e (k|p,q) + ^{f(k)-b*(k)}, (73) 

2 pq 

k+p+q=0 

(3, + 2^ 2 )// & (k) = I J] S"(k|p,q), (74) 

k+p+q=0 

where S c (k|p, q) and S 6 (k|p, q) are respectively the cross helicity and magnetic helicity transfer rates from modes p 
and q, to mode k. They are defined as 

S e (k[p,q) = S c (k|p|q) + r(k|q|p), (75) 
S fe (k|p,q) = S*(k|p|q) + S fc (k|q|p), (76) 

where 5 f (k|p|q) and S*(k|p|q) are the transfer rates from mode p to mode k, with the mode q acting as a mediator. 
Hence 

S c (k|p|q) = 3{-[u(q)-k][u(p)-b(k) + b(p)-u(k)] + [b(q)-k][u(p)-u(k) + b(p)-b(k)]}, (77) 
S*(k|p|q) = 2K{u(q)-[b(p)xb(k)]}. (78) 

We note that 5 c (k|p|q) = -5 c (p|k|q) and S 6 (k|p|q) = -S & (p|k|q), meaning that the transfers from p-to-k and k-to-p 
are opposite, as expected from mode-to-mode transfers. 

Defining the transfer rates from the shell m to the shell n, as 

7^= J] 5f ( k IPl<l)> r «m= X 5 fe (k|p|q), (79) 

kefl„,pefl,„ kefl„,pe« fJi 
q=-(k+p) q=-(k+p) 

the helicities in shell n as 

H% = J] # e (k), H b n = J] H b (k), (80) 

kea„ keo„ 
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gives 



d t H b n 



D c + F, 

nm n n 1 



V r 

m 

y T b _ D b 

/ j nm n 



with 



K = (v + 77) J] ^// C (k), = 277 J] fc 2 // fo (k), Fj = £ & {f(k) • b*(k)} . 



(81) 
(82) 

(83) 



kea„ 



ken,, 



kea„ 



The inverse cascade of magnetic helicity has been studied numerically by lAlexakis et al.l (12006b . on the basis of a 
similar shell-to-shell formalism. 

2.2.5. Helical decomposition 



Following the approach presented bv lWaleffd d 1992b for HD turbulen ce, we introduce, in Fourier space, a base of 
polar ised helical waves defined as the eigenvectors of the curl operator ( Craval 1958 ; HerringL 1974 ; Cambon and Jacquinl 



1989), 



;k x = ±kh ± 



(84) 



Note that the helical vectors h ± (k) are defined up to an arbitrary rotation of axis k. IWaleffd (119 92) suggests taking 

h*Ck) = u 2 (k) ± iui(k) (85) 

with Ui(k) = (Zk x k)/|(zt x k)| and U2( k) — Ui(k) x k/fc , where Zk is an arbitrary vector that, in general, may depend 
on k, though it is not proportional to k. iLessinnes et al.l (I2009bl) extended this approach to MHD with the following 
line of argument. 

The Fourier modes of both fields, velocity and magnetic, are expanded on that helical base 

u(k) = M + (k)h + +iT(k)h~, (86) 
b(k) = b + (k)h + +2r(k)bT, (87) 

leading to energies and cross helicity expressions 

E"(k) = i(|M + (k)| 2 + | t r(k)| 2 ), (88) 
E\k) = i(|£ + (k)| 2 + |fc-(k)| 2 ), (89) 
H c (k) = ^(u + (k)b + *(k) + u-(k)b-*(k) + c.c). (90) 

The vorticity a> and potential vector a (with an appropriate gauge) can also be expanded on the same base 

a)(k) = £(i< + (k)h + -iT(k)lT), (91) 
a(k) = r 1 (b + (k)h + -fo-(k)h-), (92) 

leading to the kinetic and magnetic helicities 

H"(k) = ^(|« + (k)| 2 -| M -(k)| 2 ) (93) 

H b (k) = &- 1 (|fc + (k)| 2 -|fc-(k)| 2 ), (94) 

(95) 
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enstrophy and square potential 



E(*) = ifc 2 (| M + (k)| 2 + | M -(k)| 2 ) 
Mk) = ir 2 (|fo + (k)| 2 + |^(k)| 2 ). 



(96) 

(97) 
(98) 



Replacing the expressions ( !86H87b for u and b in Eqs. ( 13711381 . and projecting on to the helical base h !i (k) (s k = ±1) 
leads to the following system 



(d, + vk 2 )u\\i) 
(d, + r ] k 2 )b st (\i) 



\ Z J](s P p-s q q)g(ii s i'(p) U Xq)-bHp)b s '<(q)y + / Sl (k), 

PI S„' S 1 

k+p+q=0 



(99) 



(100) 



P.<J Sp,S, 
k+p+q=0 



where g, a function of k, p, q, s k , s p and s 9 , is defined as 



g(k,p,q, s k ,s p ,s q ) = 



1 



h' !i (k)* ■ h' 5l (k) 



(h"(k)* x h s "(p)*) • h s '(q)*. 



(101) 



Considering a single triadic interaction k+p+q = 0, it is not necessary to introduce an arbitrary unit vector z k to define 
the unit vectors Ui and 112. Indeed, there is a natural direction which is represented by the unit vector perpendicular to 
the plane of the triad: 

^ = (kx p)/|k x p| = (p x q)/|p x q| = (q x k)/|q x k| . (102) 



A second unit vector fx k — k x A/k is introduced and the helical vectors are then defined as 

h*(k) = e'** (A + is kf i k ) . 



(103) 



The angle (p^ defines the rotation around k needed to transform the base (// k , A) onto the base (ui(k), U2(k)). Since the 
base (fi k , A) depends on the triad, the angle <pk is also a function of (k, p, q). The coupling constant for this triad then 
simply reduces to 



g(k,p,q,s k ,s p ,s q ) = -e 



_ e -'(.S t <f>lL + Sp'Pp+S<,(p q ) 



-i e 



-i®kpq(.Sk<Sp<<lq) 



$k s p s q (s k sin a k + s p sin a p + s q sin a q ) , 
{s k sin a k + s p sin a p + s q sin a q ) 



(104) 
(105) 



where the phase Q>k pq (s k , s p, s q) = s k {ip k + ?r/2) + s p (<p v + n/2) + s q (<p q + n/2)), and a k , a p and a q are the triad angles 
(see Fig. |5]l. 




Figure 5: Representation of the triad formed by the wave vectors k, p and q. 



The sines are defined analytically as: 



sin = — sinap = ^ sina-o = ^ , (106) 

2pq 2kq 2k p 
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where Q = y/2 k 2 p 2 + 2 q 2 p 2 + 2 q 2 k 2 - k 4 - q 4 - p 4 . The expression (1105b shows that g depends on the shape of 
the triangle formed by the triad but not on its scale. In the ideal limit (v = r] = 0) and in the absence of external 
forcing, the triadic dynamical system obtained by neglecting all the interactions with wave vectors different from k, p 
or q is given by 



d,u Sk (k) -- 


= g(k,p,q,s k ,s p , 


Sq) 


(s p p 


- Sqq) 


(wMp) u s *(q) - 


bHp)bHq))\ 


d t u s p(p) - 


= g(k,p,q,s k ,s p , 


Sq) 


{s q q 


- s k k) 


(u s "(q)u Sk (k)- 


b s 'i(q)b s "(k)y, 


d t u s i(q) - 


= g(M, p, q, s k , s p , 


Sq) 


(s k k 


- s P p) 


(u s "(k) u s "(p) - 


b St (k)b s p(p)Y, 


d t b s *Qs) -- 


= g(Kp,q,s k ,s p , 


Sq) 


(- 


s k k) 


(uHp)bHq)- 


b s p(p)u s «((l)Y, 


d,b s *{p) -- 


= g(Kp,q,s k ,s p , 


Sq) 


(" 


SpP) 


(u s "(q)b St (k)- 


u s "(q)b St (k))\ 


d t b°*{q) -- 


= g(Kp,q,s k ,s p , 


Sq) 


(" 


Sqq) 


(u St (k)b s "(p)- 


u Sk (k)b s -(p)y. 



(107) 



This dynamical system couples six complex variables. The geometric and scale independent g factor is the same in 
all equations. The second prefactors in Eqs. j 107b only depend on the wave numbers of the triad or, more specifically, 
on the eigenvalues of the curl operator. The nature of interactions in Eqs. ( 11071 ) is obviously affected by the values 
of the parameters Sk = ±1, s p = ±1 and s q — +1 (eight possible choices). However, the structure of the system is 
unchanged if all the signs of s k , s p and s q are reversed. Therefore, there are only four different types of interaction. In 
each triad (k, p, q), Eqs. ( 11071 ) automatically conserve the total energy £"(k) + E b (k) + E"(p) + E b (p) + E u (q) + E h (q), 
the magnetic helicity H b (k) + H b (p) + H b (q), and cross helicity H c (k) + H c (p) + H c (q). This is also true for kinetic 
energy and helicity in HD. 

Such helical decomposition turns out to be useful in the derivation of shell models of turbulence, mainly because 
the kinetic and magnetic helicities are then unambiguously defined. 

3. Derivation of MHD shell models 

3.1. Principles and generic equations 
3.1.1. The HD GOY model as a first example 

Shell models have been elaborated keeping in mind the spectral representation of the Navier-Stokes equations. 
The Fourier space is divided into spherical shells, defined by 

k„ < \k\ < k n+1 (108) 

for which an illustration is given in Fig. [6] 

log£ 




loglk 



Figure 6: Illustration of the shell n (in grey ). 

The sequence of k„ is chosen to be geometric with the common ratio A. Therefore k„ = koA" and the kinetic energy 
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in a given shell n is given by 

E"(k)dk. (109) 

Next we introduce a complex quantity U„, such that \U„\ 2 /2 characterizes the kinetic energy E„ in shell n. This 
quantity U„ depends on time only and is interpreted as a typical velocity fluctuation in shell n, a kind of collective 
variable for all fluctuations u whose Fourier images belong to shell n. As turbulence is assumed to be homogeneous 
and isotropic, all directions are equivalent and so all directional information is lost. 

The simplest HD shell model consists of a system of ordinary differential equations of the form 

d t U n = Q n -vk 2 n U n + F n , ne ,N\, (110) 

which mimics, to a degree, the original Navier-Stokes equations. The term -vk\U n corresponds to the viscous dis- 
sipation of U„ in shell n. The last (complex) term F n corresponds to the forcing applied in shell n. In general, F n 
is applied only in one shell. However, nothing prevents from applying F„ in several shells in order to control, for 
example, helicity injection in addition to energy injection (see Sec. 14.1.11 . 

Without the Q„ term, Eq. ( 11 lOl i is a simple diffusive equation with each shell being independent of the others. 
The Q„ term mimics the non-linear interactions within triads. A variety of shell models can be derived depe nding on 
the choice of Q„. As an introducto ry example we choose Q„ to have the form of the so-called GOY model (IGledzer , 



1973HYamada and OhkitaniLll987l) 

Q n = ik„ [aU n - 2 U„-i + bU„- X U n+ i + cU n+x U n+2 \* (111) 

where a, b and c are real coefficients. The expression of Q„ is inspired by the Fourier form of the non-linear terms 
in Eq. ( l37l i. We keep only the transfers within the subset of triads defined by (k„^2,k„-i,k n ), {k n -\,k„,k n+ \) and 
(k n , k n+ \ , k n+ 2). 

To obtain a model of 3D turbulence the coefficients a, b and c are derived writing that kinetic energy E u and 
helicity H u are ideally conserved. In this model the latter quantities are defined as 

E U =±Yj\U„\ 2 , H U = l - Yj.-mnWn?. (112) 

In absence of forcing and for v = 0, the equations d,E u = and d,H u = take the form 

£ Q«C: + CC. = 0, J^i-lTknQnUl+CC^O (113) 

where c.c. denotes the complex conjugate. This leads to 

i k„ [aA n _i + bA n + cA n+ i] + c.c. = (114) 

^(-l) n ^[aA n _ 1+ M n + cA n+1 ] + c.c. = (115) 

where A„ = £/* j U*U* +l . With appropriate subscript changes Eqs. (111411115b are satisfied if and only if 

ak„ + i + bk„ + ck„-\ = (116) 

ak^-bkl+ck^ = 0. (117) 

Replacing k„ by k^A" leads to 

a = -c/A 3 , b = -c(A-l)/A 2 . (118) 

Time and viscosity can be renormalized respectively by c and c~' , leading to an c-independent problem. Taking A — 2 
Eq. ( 11 10b becomes 



d t U„ = ik„ 



1 1 

U„+iU n+ 2 - -U n -\U n+ \ - —U n -iU„-2 



-vklu„ + F„, ne{l,--,N], (119) 
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which is the GOY model for 3D HD turbulence. 

A GOY model can also be derived for 2D HD turbulence writing that enstrophy (defined in Eq. (1132l i) instead of 
helicity is conserved. Taking A - 2 Eq. (11 1 Ob becomes 



d,U„ = ik n 



5 1 

Un+\Un+2 - -Un-\Un+\ + TT U n -l U n -2 

o lo 



v%U n + F n , b 6 {!,•■■ ,N}. (120) 



3.1.2. General formalism ofMHD shell models 

Any shell model can be recast within a general formalism such as the one introduced in lLessinnes et alj d2009al) 
and Lessinnes This formalism has the advantage of enhancing the link to the original MHD equations and of 

clarifying the differences between the variety of shell models elaborated so far. The equations are given by 

d t V = Q(U) - Q(B) - vD(U) + F, (121) 
d,B = W(U, B) - W(B, U) - t]D(B), (122) 

where U and B are vectors in space C^, N being the total number of shells, 

U = (U U U 2 ,--- ,U N ), B = (B U B 2 ,--- ,B N ). (123) 

The coordinates U„ and B„ thus correspond to the velocity and magnetic fluctuations in shell n. 
The linear operator D is defined as 

D(X) = {k\X u k 2 2 X 2 , • • • , k 2 N X N ). (124) 

The vector F stands for forcing and has non-zero coordinates only in shells which experience forcing. The operators 
Q and W stand for the non-linear terms in Eqs. ( 13711381 ). The general expression of the n-th coordinate of Q and W 
are assumed to be of the form 



Q„(X) = k n ^a%X i X j + b%X*X j + c%XX j +d%X*T j . (125) 

N 

W„(X, Y) = k n J] a^Yj + b%jX*Yj + c^Y* + d%jX*Y*. (126) 
U=l 

As an example, in the GOY model = b%. = c^. = 0, leading to the general form g„(X) = k n 2 d®X*X*. This 
choice is arbitrary. 

To determine the remaining coefficients some additional criteria have to be applied. These are: 

1 . the number of variables per shell, 

2. the number of shells interacting with a given shell «, 

3. the locality of the interactions between the shell n and the other shells, 

4. the conservation laws, 

5. the symmetries coming from Eqs. (Q]|2]i. 

In most MHD shell models, the variables in each shell n are U„ and B„. However, in helical models such as the 
one presented in Sec. 13.4. 21 the number of variables is doubled to t/+, U^, B^ and B~. 

In shell model terminology we distinguish the first-neighbor models from two-first-neighbor models, and the local 
from non-local models, depending on the interacting triads. 

• In Ll-models (local, two feet in the same shell, the third foot in a neighboring shell), the shell n is involved in 
triads (« - 1, «, n) and («, « + \,n + 1), or (« - 1, n - 1, «) and («, «, « + 1), or all the four. 

• In L2-models (local, three feet in three neighboring shells), the shell n is involved in the triads (n - 2, n — 1 , n), 
(n - l,n,n+ 1) and («, n + l,n + 2). Clearly the GOY model introduced above is a L2-model. 
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• In Nl -models (non-local, two feet in the same shell, the third foot in an inner shell), the shell n is involved in 
triads (n — m, n, ri) and (n, n + m, n + m). 

• In N2-models (non-local, two feet in two neighbor shells, the third foot in an inner shell), the shell n is involved 
in triads (n — m — 1 , n — 1 , n), (n — m,n,n+ 1 ) and (n,n + m,n + m + 1 ). 
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Figure 7: Interactions between three modes belonging to shells n, p and q. The graph is given in the map (p,q). The labels are defined in the text 
and refer to different types of shell models. The grey intensity is proportional to the probability of interaction for A = A g , the darker the higher the 
probability. White corresponds to zero probability. Adapted from Plunian and Stepanov (20071). 

This classification is illustrated in Fig. Q where the interactions between three modes belonging to shells n,p 
and q are reported in the map (p, q). The possible interacting triads are dictated by the geometry. As mentioned in 
Sec. 12.2.11 a triad is defined by three wave vectors k, p, q such that k + p + q = and therefore define a triangle as 
in Fig. |5] Here each vector belongs to a shell. For a common ratio A — 2, triads of the form (n, n,n + m) with m > 1 
are therefore impossible, as it would require a triangle with two identical sides and a third side larger than the sum of 
the two others. Of course the possible interactions depend on A. However, only the interactions corresponding to LI, 
L2, Nl and N2 are kept in shell models. The grey squares represent the relative probability of interactions between 
three shells calculated for A = A g , where A g = (1 + y5)/2 is the golden number. In this case we see in Fig.|7]that 
using a shell model four possible interactions are ignored («, «, ri), (n, n—2,n — 2), (n, n, n + 2), («, n + 2, ri). However, 
compared to the case A — 2 in which (n, n, ri) is also ignored and («, n — 2,n — 2), («, n,n + 2), (n, n + 2, ri) do not exist, 
the results do not change significantly. 

The kinetic and magnetic energies E u and E B , and the cross helicity H c are defined by 

E u = -U 2 , E B = -B 2 , H C =U-B, (127) 
2 2 

with the scalar product 

1 N 

X-Y=-Y J (XiY* + YiX*). (128) 
Z i=i 

Note that E u is homogeneous to the square of a velocity. Therefore an energy spectrum E„ = ^ U 2 obeying to a 
scaling law E% oc k" would give an energy density scaling law E u (k) oc fc^ . Typically the -5/3 Kolmogorov scaling 
law for E"(lc) corresponds to E„ oc k„ 3 '. 
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The conservation of total energy and cross helicity leads to 



[Q(U) - Q(B)] ■ U + [W(U,B)-W(B,U)]-B = 0, (129) 
[Q(U) - Q(B)] B + [W(U,B)-W(B,U)]-U = 0. (130) 

These equations must be satisfied for any vector U and B. In particular they must be satisfied when U and B are 
exchanged. This shows that both relations conservation of total energy and cross helicity are equivalent, in agreement 
with Sec. l2~L2l 

Eqs. d!29H130b must also be satisfied for B = (or U = 0) implying that 

Q(X)-X = 0, (131) 

for any variable X. 

Since a curl operator is included in the definition of the potential vector, the latter is not trivially defined in shell 
models. Therefore such a general framework fails to provide a general equation for the conservation of magnetic 
helicity in MHD or kinetic helicity in HD turbulence. Actually the definition of the latter quantities depend on the 
type of shell model used, helical or non helical. It is postponed to Sec. l3.2l[3.3l and l3~4l The enstrophy S and squared 
magnetic potential A are more easy to define as, contrary to helicity, they are always positive in any shell 

2 = \ X ^ 2|t/ » |2 ' A = \ X fc " 2|B » |2 - (132) 

n n 

We note that the necessity of having the same structure LI, L2, Nl or N2, for Q and W follows from the conser- 
vation laws. 



At this stage we have all the information necessary to start elaborating a shell model. It is a matter of choosing 
between the several possibilities mentioned earlier. However, we note that contrary to Q, the operator W is not 
uniquely defined. Indeed W(X, Y) can be replaced by any operator of the form 



W(X,Y) 



V -W(X, Y) + — !— W(Y,X) 



(133) 



w — 1 w — 1 

where w is a scalar quantity, without changing Eq. ( 1122l i. or the conservation laws. Indeed it is easy to show that 

W(X, Y) - W(Y, X) = W(X, Y) - W(Y, X). (1 34) 

In other words W(X, Y) corresponds to a combination of (x ■ V)y and (y ■ V)x. Though this does not change the shell 
model results in terms of U and B, it may change their analysis (post-processing) in terms of energy transfer and 

shell-to-shell exchange. 

Following lLessinnes et al.l d2009al) we arbitrary make W(X, Y) correspond to -(x ■ V)y. From Eq. ([8]) this implies 

that 



W(X,Y)-Y = 0, 



and uniquely determines the value of w in Eq. ( 1133b . 

So at this stage Eqs. (1 1 2 1 H 1 22b can be rewritten in the form 



d,V = Q(U) - Q(B) - vD(U) + F, 
d t B = W(U, B) - W(B, U) - /7D(B), 



(135) 



(136) 
(137) 
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with 



Q(X)-X = 0, (138) 

W(X,Y)-Y = 0, (139) 

Q(X)- Y + W(X,Y)-X = 0. (140) 

W(X,X) = Q(X), (141) 



which again is present in the original Navier-Stokes and induction equations. To show that Eq. (11411) is satisfied, we 



start by replacing Y by X + Y in Eq. ( 1139b . Assuming that W is a bilinear operator, 

W(X, X) ■ X + W(X, X) ■ Y + W(X, Y) ■ X + W(X, Y) ■ Y = (142) 
which, from Eq. ( 11391 ) again, simplifies to 

W(X, X) ■ Y + W(X, Y) ■ X = 0. (143) 

From Eq. ( [T40l W(X, Y) ■ X can be replaced by -Q(X) ■ Y in Eq. (fl42l leading to 

(W(X, X) - Q(X)) ■ Y = 0. (144) 

As Eq. (1144b must be satisfied for any Y, it implies Eq. (11411) . 
Finally Eqs. J 12 1H 122b can be rewritten in the form 

d t V = W(U,U)-W(B,B)-vD(U) + F, (145) 

d t B = W(U,B)- W(B,U) -rjTHB), (146) 

with 

W(X,Y)-Y = (147) 

being the only condition required to be satisfied so that both total energy and cross helicity are conserved along with 
the symmetries of the non-linear operators in the original Navier-Stokes and induction equations. In particular the 
Eqs. ( 11381 ) and ( 1140b are automatically satisfied, using again the bilinear property of W to show Eq. ( 11401 ). 

Note that changing U to KXJ, B to KB and W to K~ l W, where K is scalar quantity, does not change the system of 
Eqs. (I145H146I) . 

Similar to W, W takes the general form 

N 

W n (X, Y) = J] aZjXiYj + bZjX*Yj + cf^X/Y* + d%X*Y*. (148) 
U=l 

As the velocity u and induction b are divergence-free, the Navier-Stokes and induction equations satisfy Liou- 
ville's theorem in the ideal limit v = r\ — 0. For shell models the Liouville theorem gives 

y _d_ IdVA _d_ IdBA _d_ tdu:\ + J_(d_K\ 

t-i dU„ \ dt / 3B n \ dt ) dm \ dt ) 8B* \ dt I 



Finally, the magnetic helicity which has not been considered so far will give rise to an important additional 
constraint on shell models designed for 3D MHD turbulence. 
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3.1.3. Energy Flux 

Starting from Eqs. (114511146b we obtain the kinetic and magnetic energy equations in any shell i 

d t Uf/2 = Tf u + T? b -vIcfuf+F-Vi, 

BB , T^BU 



d,Bf/2 = Tf" + Tf u -jjkjBj, 



where 



Tf u = W(U, U) ■ Vi, Tf B = -W(B, B) • U;, 



W(U,B)-B,-, Tf 



-W(B,U)-B,-. 



(150) 
(151) 



(152) 
(153) 



and X, = (0, ■ ■ ■ ,0, X,, 0, ■ ■ • ,0). The quantities Tf Y are interpreted as the energy rate flowing from all shells of the 
Y-field into the ;-th shell of the X-field. They are illustrated in Fig. [8] 



N 





1 ■■■ i-1 i i+1 ... N 1 j N 

Figure 8: Illustrations of the energy rate flowing into the i'-th shell of X, coming from either all shells of Y (left) or only one shell j of Y (right). 

Eqs. (fT43l and imply respectively £ Tf B + £ Tf u = and f TV U = f Tf B = 0. 

i=l ' i=l ' i=] ' i=l 

Going one step further we introduce the quantity Tf?, the shell-to-shell energy exchange rate from shell j of the 
Y-field into shell ; of the X-field, defined as 



T™ = W(U, U,-) ■ Vi, T" B = -W(B,B,-) ■ V u 
Tf B = W(U, By) • B,-, T BU = -W(B, Vj) ■ B h 



implying that 



nXY 



7=1 



A')' 



that is illustrated in Fig. [8] Then the energy exchange rate from Yj to X, must be opposite to that from X,- to Yf 



'-pXY '-pYX 

U ~ ji ' 



(154) 
(155) 

(156) 
(157) 



This implies the conservation of total energy (1 1 29b . the relation (1131b for X = U, and Eq. (1135b . The notation Tf? is 
chosen in agreement with T„ n given in Eq. d64l l. 

In order to calculate energy fluxes we introduce the following vectors 



X< = (Xi, ■ ■ ■ ,X n , 0, • • ■ , 0), 
X> = (0,--- ,0,X n+u --- ,X N ), 

with X = X,^ + X^. The corresponding energies are defined in the same manner as in Eq. (1127b . 

Ef = l -(Xf, Ef = i(X„>) 2 . 



(158) 
(159) 
(160) 



(161) 
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Then from Eqs. ( I150I15U . 



d t E? = ^(r,^ + r,™)-vD(U)-U,f +F-U< 

;=i 

d,Ef = J(r, M + rr)-^D(B)-B„<, 



(162) 
(163) 



with D(X) ■ X^ = 2 kfXf . The quantity £ Tf Y defines the flux from Y (in all shells) to X< It is denoted by 



We also define 



(164) 



!=1 



nf< - Z Z T " 



zz^ 

1=1 J=l 

ZZ 7 ^ n *-ZZ^ 



i=l 7=n+l 



xy 



(165) 



i=n+l 7=1 



(=n+l j=n+l 



which is the shell model counterpart of Eqs. d69ll72l i. In the flux notation given in Eqs. dl 64111 65b the subscript n has 



>'„ 



been dropped for convenience, e.g. flL must be understood as II ". We have 



u y x< = n£+n£, nf =n£+n£ 

nT = n£ +ni:. 



n£ + n£, 



(166) 



Using Eq. (I157l i we can show that 



n z< - 
ii x < - 


= n£ = o, 








= -n£, 


nf< = 


-n x< 




- -n£, 


n£ = 





(167) 




Figure 9: Illustration of some energy fluxes. The others can be found with the help of Eqs. U67I . 

These equations, illustrated in Fig. [9] imply the following expressions for the energy equations 

d t E^ = n£+n£+n£-vD(U)-i£+F-i£, 



d,Ef = n£ + n£ +n u B : -77D(B) b<, 



(168) 
(169) 
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where 



n£ = w(u,i£)-u<, n£ = -w(b,b<) u^, 



U B V : = -W(B,B>)-U<, 
n£ = -W(B,U„<)-B<, 



Il£ = W(U,B>)-B<, 
n£=-W(B,U„>)-B< 



The evolution of kinetic and magnetic energy in shell i has the form 



N 1 W 1 

J ( t/, 2 /2 = J] -S uu (i\j, k)+Y, ^S UB (i\j, k) - vk) Uf +F ■ U jt 

d,Bj/2 = j^ l S BB (i\j,k) + j] l -S BU (i\j,k)-^Bj, 



(170) 

(171) 
(172) 



j,k=i 



j,k=i 



with 



S uu k) = W(U,, U,) ■ U< +W(U,-, U*) ■ u jf 
-WCBfaB^-U, -WCBj.BjO-U,-, 
W(Ut,B;)-Bj +W(ty,B*)-Bi, 



S ra (iL/,*) 

s M (iL/,*) 



S au (i|j, it) = -W(B t , U,) ■ B, -W(B ; , U*) • B,-. 

Each S XY (i\j, k) term represents the transfer rate from the modes j and k of field Y into the mode i of field X. 
We also define the following quantities 



(173) 



S uu (i\j\k) = 
S UB (i\j\k) = 
S BB (i\j\k) = 
S BU (i\j\k) = 



W(V k ,Vj)-Vi, 
W(B k ,U,)-B jt 



(174) 



where S XY (i\j\k) is the mode-to-mode energy transfer rate from the mode j of field Y to the mode i of field X, with the 
mode k acting as a mediator. Within one triad (i, j, k) the related interactions are illustrated in Fig. [10] The following 
relation applies 

rf = JV r 01#). (175) 

k 

MHD fluxes were intro duced in lStepanov and Plunianld2006l) and later corrected in lPlunian and Stepanovl (120071) and 



Lessinnes et al. (2009a). 



3.2. Local models 

3.2.1. LI -models (local, first-neighbor) 

In Appendix A| we give the form of all possible Ll-models ob eying to the general req uirements given by Eqs. ( 1147t 
and (1148l l. They were al ready introduced in th e seminal paper bv lGloaguen et al.l (11985b. Two Ll-m odels have been 
studied in detail, one by Gloaguen et al.l d 1985b and subsequent authors, the other by Biskampl l 1994 ). 



The model investigated bv lGloaguen et al.l (119851) corresponds to 

W n (X f Y) = k n [diX^Y^ - AX„Y„ +l ) + C 2 (X„Y„^ - AX n+1 Y n+l )] , (176) 
where C\ and C2 are real parameters, and X and Y are real variables. The set of equations for U and B are as 
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Figure 10: Illustration of the mode-to-mode energy transfers. The grey triangles correspond to the interacting triads. 



follows 



d,U„ - C\ (k n U 2 _ x - k n+ iU n U n+ \ - k n B\_ x + k n+l B„B„ +l s j 

+ C 2 (k„Un-i U n - k n+1 U 2 n+l - k„B„^B„ + k n+1 B 2 n+1 ) - vk*U n + F H , 
d t B n = C\k n+ \ (U n+ \B n - U„B n+ i) + C^kn (JJ n B n -\ - U n -\B n ) - r]k\B n . 



(177) 
(178) 



In addition to total energy and cross helicity, provided C\ +A q Co = 0, this model has another conserved quantity, 
2 k\B n . However, as it is not quadratic it cannot be considered as an analog of magnetic helicity (which is 
conserved in ideal 3D MHD) or as a squared magnet i c pote ntial (which is conserved in ideal 2D M HD). So it 
has no real meaning. For B = 0, the Bell and Nelkin (I 1978b model is recovered, which in turn gives Obukhov 
dl97ll) model if C x = 0, and lDesnianskii and Novikovld 1974 model if C 2 = 0. 



In the dissipationless limit (v = r\ — 0) and for an infinite number of shells, the system of Eqs. (1177H178t has 
the Kolmogorov stationary solution 

B„ ~ U„ ~ k~ U3 (179) 



for any value of the ratio c = C1/C2. However. iGrappin et al.l (119861) found three kinds of attractors depending 
on the value of c. For c > 2 the system tends to a supercorrelated state, for which the attractor is reduced to a 
fixed point B = ±U, characterized by a total absence of spectral energy transfer and very steep energy spectra. 
For 0.5 < c < 1, the attractor is a non-magnetic (B = 0) stable fixed point. Finally, for c < 0.3 the attractor has 
a high dimension, with a chaotic solution characterized by an extended inertial range and equipartition between 
kinetic and magnetic energies. Therefore only the latter range of c is of inter est for MHD turbule nce as it 
reproduces the expected chaotic behavior of real turbulence. Finally, for c = O.Ol. lGrappin et al. (1986) found a 
Lyapunov dimension for th e attractor which is consistent with the standard Kolmogorov HD turbulence, rather 
than the Kraichnan MHD ( Kraichnanl 1965 ). phenomenology . The latter scenario is lost due the absence of 
non-local interactions in the model. 



In its HD version, i.e. for a model reduced to Eq. ( 1177l i with B = O. lDombre and Gilsonld 19981) found a solution 
for U which again depends on c. A stable fixed point is found for c > 0.55 (including Desnianskii-Novikov's 
model for c — > 00), and chaotic solutions for c < 0.536 (including Obukhov's model for c = 0). The transition 
between both regimes corresponds to a succession of Hopf bifurcations. 



MHD intermittency has been studied by ICarbone Jl994alTbh for c = 0.01 and = 19. Solving the system of 
Eqs. (1177H178t expressed in terms of Elsasser variables Z* = U„ + B„, he calculated the p-th-order structure 
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functions S*(p) = (\Z*\ P ). For scales belonging to the inertial range, he found a scaling exponent ( p such that 
S„(P) K V- For a range of scales much larger than the inertial range, he found that the structure functions 
satisfy the relation S*(p ) oc S^{3Y P ^\ thus confirming that the concept of extended self-similarity applies to 
MHD. ICarbond dl994bl) also compared thes e scaling expon ents to those obtained from solar wind measure- 
ments, by the Voyager 2 satellite at 8.5 AU (IBurlagai 1 1 99 II) . As show n in F igfTTl the shell model results lie 
inside the error bars of the observed data. For completeness we note that lCarbond ( 11995b showed that the solu- 
tions of Eqs. ( 1177lll78t are sign-singular, meaning that their sign reverses cont inuously on arbi trary finer time 
scales, similarly to some signed measurements in turbulence and fast dynamos (lOtt et al.[ll992l) . 




2 4 6 8 10 12 14 16 18 20 22 



Figur e 11: Scaling exponents versus p obtained from th e Ll-model U76t (full line), and from Voyager 2 solar wind data analysis jBurlagal . 

1991) (points with error bars). Adapted fromjCarbone 1 1994^). 



The model investigated by Biskampl (1 1994 ) corresponds to 

w n (x, Y) = k n [c.ixiX-, - *KK+d + WtCi - CCi)] - 

and gives 



(180) 



d,u„ = d (k n u* n : x - k„ +1 u*„u* n+l - k n B* n \ + k„ +1 B;x + i) 

+ C 2 (KUUK - k n+ xUll x - KK-xK + k„ + iB* n 2 +x ) - vk\\J n + F n , (181) 

d,B„ = C^^B^-^B^ + Cz^^X-i-^-iS:)-^, (182) 

where U and B are complex vectors. Again, depending on the ratio c - Ci ICt, iBiskarnpl ( 1994 ) studied the 
dynamics of the solutions of Eqs. ( 118111182b for both HD and MHD turbulence. 

For HD t urbulence, chao t ic Kol mogorov solutions were found only for |c| «: 1 . This confirms the argument 
made by iGloaguen et alj ( 11985b that \c\ «: 1 is more representative of incompressible turbulence, since the 
C\ -terms referring only to flat triads make a negligible contribution to the non-linear interactions. 

For MHD turbulence Biskampl ( 1994 ) calculated the structure functions scaling exponents of Z*, U and B 
for c = -1.33 and c = 10~ 2 . He found the same scaling exponents for all variables (within the error bars), 
but the results depend on c. For c = -1.33 he found stronger multifractal behavior (stronger deviation from 
the Kolmogorov scaling exponent £ p = p/3) while f or c = 10 ~ 2 he f ound scaling exponents that were more 
compatible with experimental observations and DNS. Biskampl ( 1994 ) also studied the effect of an externally 
applied magnetic field. 
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It is remarkable that th e peak in popularity of LI -models reached at the beginning of the n ineties was mainly due to 
the Gloaguen et al. ( 1985 ) model which, in contrast to the HD models of Obukhov ( 197 lb and Desnianskii and Novikov 
( 1974 ). demonstrated not only chaotic behavior but also intermittency. This lead Brandenburg! ( 19921) to generalize 
this model to MHD Boussinesq convection by adding an equation for tempera ture fluctuation 6. For B = he f ound 
(for a range of parameter c) the scaling laws E u {k) ~ A:"/ 5 and Eg(k) ~ k^ 1 ^ 5 ( Obukhov , 1959tlBolgianol T959I). For 
B + he found Kolmogorov spectra for the three fields Eg(k) ~ E u (k) ~ E^ k ) ~ k~ 5 ^ 3 . iGeertsema and Achterberg 



01992b introduced a vectorial three-component version of the iGloaguen et al.l d 1985b model in order to evaluate the 
turbulent stress tensor in a differentially rotating disk. This model will be discussed in Sec. 14.3.21 



A generic problem encountered with the previous model s is that there is an insufficient num ber of quadratic 
invariants. In the HD models introduced by lObukhovl ( 1197 lb or iDesnianskii and Novikovl d 1974b only the kinetic 
energy is conserved bringing th e number of quadratic inv ariants to one, in stead of two in real HD turbulence. In 
the MHD models introduced by Gloaguen et al. ( 1985 ) or Biskampl ( 1994 ) only the total energy and cross helicity 
are conserved, bringing the number of quadratic invariants to two, instead o f three in real MH D turbulence. Such a 
problem is due to a too simplistic expression of W„(X, Y). For example the Biskampl ( 1994 ) model corresponds to 
only the four first terms in Eq. ( lA.2t among the twenty possible terms (assuming A io = 0). However as shown in l3.4.1l 
one additional quadratic invariant can be obtained taking more terms in Eq. dA.21 ). Such a model is called helical and 
will be detailed in Sec. 13.4.11 It is called helical for at least two reasons. First it allows to have the kinetic helicity 
in HD or magnetic helicity in MHD as additional quadratic invariant. Second such a model can be interpreted in a 
framework of helical mode decomposition. 



3.2.2. L2-models (local, two-first-neighbor) 

Another way to introduce an additional quadratic invariant is to increase t he number of in teracting triads, with 
interactions between two-first-neighbors instead of just first-neighbors. This lead lGledzerl ( 1 1 973b to propose a "System 
of hydrodynamic type admitting two quadra tic integrals of motion" (his paper's title). In our notation it is a L2-model. 
It is remarkable that at about the same time lLorenzl d 1 972b derived exactly the same model, starting from the Navier- 
Stokes equations. Both authors imposed enstrophy conservation in addition to kinetic energy, thus the model is 
relevant to 2D-turbulence only. 



Following a different approach. lFrikl (11983b elaborated a non-local L2-model with again kinetic energy and enstro- 
phy conservation. Applying the same method iFrikl (11984b also derived the first MHD shell model, using the square of 
the magnetic potential as the third quadratic invariant. This model was thus relevant to 2D MHD turbulence. In our 
notation both models are N2-models. For this reason they will be discussed in Sec. 13. 31 when dealing with non-local 
models. 

Finally, Brandenburg et al] ( 1996b introduced the first 3D MHD L2- model using magnetic helicity as the third 
quadratic invariant (see also lBasu et alj dl998b : lFrick and Sokolofll (11998b ). 



The L2-model wh i ch ha s received most attention is undoubtedly the GOY model, named after iGledzerl dl973l) . 
Yamada and Ohkitani Jl987lR Among all possible L2-models given in | Appendix B| it corresponds to 



w„(x, y) = ik n [d(J5_ 2 Ci - ^C-iCi) + C2PC-1C2 - + W +1 Ci - ^C +2 Ci)] ■ d83) 

The HD version of this model became popular because of its relevance to real turbulence in terms of high-order 
structure function s. In Sec. 13.3.21 we will also discuss another L2-model, called the Sabra modeH, introduced by 
L'vov et al.l d 1998b as an "improved" version of the GOY model. It is an improvement in so far as some spurious 
correlations existing between different shells in the GOY model, are suppressed with the Sabra model. Such spurious 
correlations do not exist in the MHD version of the GOY model, both models giving the same results. 



4 Unfairly ignoring Lorenz's contribution. 

5 The name Sabra model was introduced bv lL'vov et all d 19981) in the context of shell models, presumably as an insider joke: Sabra denotes a 
Jewish person born in Israel, while Goy denotes a non- Jewish person. 
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To derive the GOY model for HD turbulence, we write W„(X, X) in the form 
W„(X, X) = -ik„A(C 3 + AC 2 ) 



A n+l A «+2 _ "T A «-1 A «+1 



1 - R 

-p- A «-l A n-2 



(184) 



where s = (C3 - AC\)/{Cj, + /IC2). With an approp riate renormal i sation of U and W in Eq. ( 1145b . the term 
-MCt, + AC2) can be taken equal to unity, leading to dBiferale et all 1 19950 



+2 



1 -e , , 



vklU n +F„. 



(185) 



Kadanoff et al.l (1 19951) showed that, in addition to kinetic energy, this system has the following quadratic invari- 
ant 

H = ^|[/„| 2 ( e -ir". (186) 



Depending on s, different physical interpretations can be given to H dFrick et aU 1 19951) : 
- For e = 1 - AT 1 , 



// = 2(-l)"fc„|[/„| : 



(187) 



which is analogous to kinetic helicity (IKadanoff et al.lll995l) . This corresp onds to the GOY model for 3D 
HD t urbulence introduced in Sec. 13.1.11 It has been studied numer ic ally (jjensenL 1991 ; Pisarenko et al 



19931) and analytically using a closure model (IBenzi et al.L Il993al) . iFrick et al. Nl995h showed that the 
GOY model displays the same intermittency as 3D turbulence provided s « 1 - AT 1 . 

- For s — 1, all cascades are impossible. 

- For s = 1 + A' 2 , 

H = Y J k L n \U n \ i (188) 

which is analogous to enstrophy dLorenz , ll972tlGledzeil 1 1 9731) . The two stable solutions are U„ oc k n 
and U„ oc k~ l , leading to spectral properties analogous to 2D turbulence. The energy spectrum density is 
indeed made of two power laws depending on whether k is smaller or larger than kf , the forcing wave- 
number. 

* For k < kf, E(k) oc k~ 5 ^ with an inverse energy cascade. 

* For k > kf, E(k) oc k~ 3 with a direct cascade of enstrophy. 

- For s — 2 the only quadratic invariant is the kinetic energy, leading to a Kolmogorov fixed point. 

- For e = 1 + A, 

H = Y J KW, (189) 

which is the dimensi onal equivalent of the "action ", a hidden integral of motion in 3D turbulence written 
in Clebsch variables ( lYakhot and ZakharovHl993h . 



For < s < 1 and A — 2, a numerical study of the transition to chaos was performed bv lBiferale et al.l i 1995b . 
They found a Kolmogorov stable fixed point solution for e < 0.3843. For s = 0.3843 the solution becomes 
unstable via a Hopf bifurcation. For larger values of e, the system evolves towards a chaotic state following 
a Ruelle-Takkens scenario. For e > 0.3953 the dynamics are intermittent with a positive Lyapunov exponent. 
This regime is characterized by a strange attractor remaining close to the Kolmogorov unstable fixed point. 



• In MHD, taking W„(X, Y) defined by Eq. dl83l . and provided 

-1 +(1 -£) + (! -e) 2 



Ci 



1 + (1 - s) + (1 -e) 2 



■c 2 , 



c 3 



, l+(l- e )-(l-e) 2 „ 
1 +(1 -e) + (l -e) 2 



(190) 
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we can show that Eq. ( 11461 ) has an additional quadratic invariant 

/ = J]( e -l)"|B n | 2 . 



(191) 



Using the same renormalisation applied to the HD case, and thus omitting the term -A(C-} + AC2) in front of 
W„(X, Y), the latter becomes 

w„(x,Y) = 1 J^[x; + x + 2 + K + 2y: +1 + { -^^(K + X + 2-K + 2y: +1 ) 



B 1 — B 

~ ^(^71-l^n+l + X n+l Y n-0 + ^(2 - E ^ X >'-1 Y «+1 ~ X n+l Y n-0 

1 - s 1 

~~ ^2 ( X n-2 Y n-l + X n-\ Y n-2) + ^2(2 _ E ^ X "-2 Y n-\ ~ X n-\ Y n-2> ] 

It leads to the system dFrick and Sokolofn,ll998tlAntonov and FrickLl2000h : 



(192) 



d,u n = ik n {{u: +l u* n+2 - b: +1 b: +2 ) - ~ etc, k +1 - bub: +1 ) - { l^-(u* n _ 2 u* n _ x - B* n _ 2 Bui\-<u n + f„, 

(l-e) /r „ - - ... . 1 



d,B n = 



2-s 



(1 - s)\U* n+x B* n+2 - B* n+l U* n+2 ) + ^1(^,2^, - B;.!^!) + ^(f/:_ 2 5;_, - B^iOj-^ 



Similar to the quantity H in HD, different physical interpretations can be given to I depending on e : 
- For s = l- 



(195) 



which is analog ous to magnetic helicity (Brandenburg et al. Il996h and, contrary to the choice of s = 
1 + A~ l made bv lBiskarnpl (119941) . leads to an unsigned quantity as expected for magnetic helicity. This 



corresponds to the GOY model for 3D MHD turbulence. 
- For s = 1 + A' 2 , 



(196) 



which is analogous to the square of the potential vector. This corresponds to the GOY model for 2D MHD 
turbulence. 



- For b = 2, 



= ^ \Bj 



(197) 



which is twice the magnetic energy. As total energy is conserved this corresponds to separate conservation 
of both kinetic and magnetic energy, which has no obvious physical meaning. 



Antonov and Frickl d2000t) solved Eqs. (|193H194t for s e [-10, 10], A = 2, v = rj = 10~ 9 , and for a stationary 
forcing applied at shell n = 0. They found that in contrast to the HD case there are no stable solutions and the 
behavior is always stochastic over the whole range of b. Kinetic and magnetic spectra are shown in Fig.[T2lfor 
different values of e. For s < 1 .01 (left column) a small-scale dynamo action occurs with near equipartition be- 
tween kinetic and magnetic energies and approximately Kolmogorov spectral slopes. For b > 1.1, the magnetic 
energy at large scales is depleted (right column). On increasing e the peak of the magnetic spectrum moves to 
smaller scales until it reaches the dissipation scale for e as 1.7. For s 6 [ 1 .7, 2] the magnetic spectrum is much 
lower than the kinetic spectrum, but the dynamo still occurs. For s > 2 the small-scale dynamo is lost. However, 
the magnetic energy at the largest scales slowly grows. To understand why it is so, we can write I in the form 
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7 = 2 kf t \B„\ 2 . Then taking s > 2 corresponds with having d > 0, 1 becoming a magnetic analog of generalized 
enstrophy. Similar to 2D HD turbulence, energy can only be transferred towards large scales (inverse cascade) 
while I is transfered towards small scales (direct cascade). 
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Figure 12: Kinetic (white) and magnetic (black) shell energy vs the shell number n for different values of the parameter e. Left(up-down) 
e = -5,-2,-0.5,0,0.5,0.99, 1.01. Right (up-down) e = 1.1, 1.3, 1.5, 1.7, 1.99,3,5. From lAntonov and Fricld f2000l) . 



Frick and Sokolofn dl998l) investigated the model d 1 93H 1 94-b in detail for s — 1 - A and e = 1 + X , corre- 



sponding respectively to 3D and 2D MHD turbulence, with I being respectively the magnetic helicity and the 
square of the magnetic potential. In free decaying turbulence (without forcing), they found dynamo action in the 
3D case, and magnetic decay in the 2D case (as expected from antidynamo theorems). In 3D (FigfOl-left) the 
magnetic energy reaches equipartition after 20 turn-over times. In the kinematic approximation, W(B, B) = 0, 
the growth of magnetic energy is unbounded as expected from kinematic dynamo action. In 2D (Fig[T3lright) 
the magnetic energy for both non-linear and kinematic cases grows up to a level of about 1/100 of kinetic energy 
and slowly decays on a dissipation time scale. 

In free decaying turbulence. iFrick and Sokoloffld 19981) also te sted the analogy between magnetic energy in 2D 
MHD turbulence and temperature gradients (IZeldovich[[l956l) . They considered a GOY model for temperature 
fluctuations 9„, with 2 8% as an additional ideal invariant. For both 2D and 3D turbulence they found that the 
thermal energy decays smoothly, while the temperature gradients exhibit temporal growth followed by decay 
similar to the magnetic energy in 2D MHD turbulence. 

Forced MHD turbulence has been tested against the concept of extended self-similarity ( Basu et all 1998b. 
while the scaling properties and long-time behaviour of the solutions have been studied in detail (IFrick and Sokolofl 
1998b . In particular, it was shown that the model given by Eqs. ( 1193H194t is rather sensitive to the dynamics of 
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Figure 13: Free decaying 3D (left) and 2D (right) MHD turbulence: kinetic E u and magnetic E B energy versus time. Thick lines corresponds to 
kinematic simulations imposing W(B,B) = 0. FromlFrick and Sokoloff 1 1 998ft . 



the applied forcing. For a constant forcing lGiuliani and Carbond (119981) found that after some time a supercor- 
relatio n state B = ±U app ears which is spurious and has no real physical interpretation. Instead of applying a 
forcing lFrick et al.l Q200CM) imposed the modulus of the complex velocity at shell n — letting its phase evolving 
freely. They found different time stages with both low or high correlations be tween U and B. Using forci ng 
with a random complex phase suppresses the problem of supercorrelation state (IStepanov and Plunian , l2006h . 



3.3. Non-local models 



Some time after the emergence of shell models in Moscow, mainly in the w ake of|O bukhov (11971ft . another 
approach, the so-called hierarchical approach, was developed in Perm (Russia) by ZiminT i 1981 ). The idea was to 
model HD turbulence as a network of vortices with a double repartition in both physical and Fourier space. Under 
the assumptions of homogeneity and isotropy, it is possible to calculate the Reynolds tensor for each interacting triad 
of vortices, and derive a system for the intensity of each vortex still dependin g on space a nd scale. Averaging in 
space over all vortices having the same scale leads to a non-local shell model. Frikl ( 1983 ) found an original way 
of enabling such a model to conserve two ideal invariants. Considering kinetic energy and enstrophy conservation, 
he elaborated a non-local shell model of 2D-turbulence. Applyi ng the same method to MHD, with conservation of 
total energy, cross helicity and the square of magnetic potential, iFrikl d 19841) developed a non-local shell mo del for 
2D MHD turbulence. Both shell m odels are N2 models. Completing the work started about 15 years earlier by Ziminl 
( 1981 ), Zimin and Hussainl ( 1995 ) derived an N 1 -model of 3D HD turbulence, which satisfies both kinetic energ y and 
helicit y requirements (though helicity is not mentioned in thei r work). The hierar c hical approach followed by iFrik 
( 1983 ) is detailed in Sec. 13.3. ll along with the model derived by Zimin and Hussain ( 1995b . 

Another way to derive a non-local shell model is to begin directly from the shell model structure given by 
Eqs. J 148b . including all possible non-local interactions. The general shape of W for Nl and N2 models is given 
respectively in Appendix C and Appendix D One example of an N2-model is presented in Sec. 13.3.21 Such direct 
derivation does not provide unique definitions for the non-local coefficients. One free parameter, directly related to the 
strength of the non-local interactions, remains. This free parameter can be estimated either from phenomenological 
argument or using the hierarchical approach. 



3.3.1. N2-model derived using a hierarchical approach 

Following IFrikl (119831) . we consider a network of parallel 2D-vortices, depending on the horizontal coordinates 
only, with velocities perpendicular to the third direction e z . Each vortex is denoted by two integers n and N. The 
first integer n indicates the shell to which the vortex wave number belongs. As with shell models we consider shells 
obeying a geometric sequence, here with A — 2 as their common ratio. 

As two different vortices may belong to the same shell in Fourier space, a second integer N is necessary to 
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differentiate the vortices in real spaced Then we can write the velocity field as 

u(f, x,y) = 2^ U nN (t) u„ N (r - r nN ), (198) 

where U„n(t) is the amplitude of the vortex nN, r„m is the position of the vortex center and u„n is defined by its Fourier 
coefficients 

<W(k) = i 2 ' " k ^ C - : e -' kr " w , for k„<\k\<k„ +l , with k„=n2", (199) 
u„jv(k) = outside the shell. 




Figure 14: Left: illustration of logarithmic shells in Fourier space, where the grey annulus corresponding to shell n. Right: functions i«„ai(a') and 
0>nN( s ) for n = 0. 

Taking the inverse Fourier transform of fi„#, and calculating the vorticity intensity io„n — (V x vl„n) • e?, we find 

u„ N (r - r„ N ) = (e z x s)(3tt)- 1/2 (/ (2j) - / 0)) /a 2 , (200) 
cj nN (r-r nN ) = 2-"(7r/3) 1/2 (27,(2^-7i(^))/^ (201) 

where s = 7r2"(r - r^), s = |s|, and 7o(5) and Ji(s) are the zero and first order Bessel functions. Fig.[14]shows both 
functions, U h n (defined such that = (e z X s)U„n) and a>„^, along with the corresponding shell n (in grey) in which 
fl n jv(k) is non zero. It can be shown that the density of vortices (number of vortices divided by the shell surface) 
increases with n as 2 2 ". 

All functions u„/v(k) and u,„M(k) are orthogonal provided n + m (functions of different scale do not overlap in 
Fourier space). This implies that u„n and u,„m are also orthogonal. In contrast, u„n and u„m are not necessarily 
orthogonal, which will have several consequences as discussed below. 

By replacing the expression for u given by Eqs. d 1 98b and ( 1200b in the Navier-Stokes equations, we find that 
amplitude U„n has the form 

d t U nN = 2^ 2 RnNpPqQUppUqQ - vkfontf- (202) 
P.P q,Q 

where non-linear interactions are given by 

RnNpPqQ = J u„ N ■ [(u pP ■ V)u ?e ] dr. (203) 
From Eq. (1203b we obtain the exact relations 

RqQpPnN = -RnNpPqQ, RpPqQnN = —RnNqQpP, RpPnNqQ = —RqQnNpP, (204) 



6 For convenience, and in this section only, N denotes the second integer and not the maximum number of shells. 
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which imply energy conservation, where energy is defined as 2 Uf lN /2. We note that this definition of energy is not 

n,N 

exact due to the fact that the u,,^ functions are not orthogonal with re spect to N (Frikl 1983b . Such a hierarchical 
model has been used to study enstrophy intermittency in 2D turbulence ( Aurell et all 1994b ). 
Now the aim is to reduce Eq. ( 12021 i to a shell model of the form 

d,U„ = J] T npq U p U q - vklU n , (205) 



p.q 



where the coefficients T„ pq need to be defined. 



We introduce the vectors R„ pq , which we assume, as in Eq. ( 1204b . satisfy 

Rqpn — Rnpq-> Rpqn ~ ~~ R-nqp^ Rpnq ~ ~ Rqnp- (206) 

Their moduli are defined to be the root-mean-square value of R„N P pqQ, calculated for all possible positions of any 
three interacting vortices nN, pP and qQ, 

|R; W | 2 = J J Rl NpPqQ dr pP dr qQ . (207) 
Frikl (119831) introduced the correlation between all triads of vortices belonging to shells n, p and q 



I J R nNpPqQRnNqQppdr p pdr qQ 

cos 6 npq = , (208) 

l^-npql \^nqp\ 

and found that 6 npq + 9 qnp + 9 pqn — 2n, and 8„ qp + 8 qp „ + pnq = 2n. Together with relation (1206b this shows that 
the six vectors R„ pq , R nqp , R pnq , R pq „,R qpn ,R q „ p are coplanar. The angles uniquely define the mutual positions of 
these vectors, but not their absolute positions as the set of vectors can be rotated by any angle. This degree of freedom 
corresponds to some arbitrary coefficient that will be set to unity after renormalisation of the equations. 



Now we define the coefficients T npq as 



T npq = (e • R npq + e • R nqp ), (209) 



where e is a unit vector. To determine the direction of e, iFrikl (119831) considered enstrophy conservation, relevant to 
2D turbulence. From the enstrophy equation he introduced the non-linear interactions 



= J w„iv(u 



V)aj qQ dr, (210) 



OnNpPqQ - J <^>nN\"pP V )U> q Q 

and the six vectors S npq , S nqp , S pnq , S pqn , S qpn , S qnp . The S-vectors are found to be coplanar with the R-vectors and the 
T n pq coefficients are given by 

Tnpq — (C ' Snpq + e ■ S n qp)- (211) 

There is, however, only one possible choice for e such that both definitions d209t and (121 U give the same value for 
T npq . For this direction of e, both energy and enstrophy are ideally conserved. 
The corresponding shell model has the form 

U n +mU„+m+l] + F„. (212) 

m>0 

The Topq coefficients were directly estimated by Frikl ( 1983 ) for X = 2, their numerical values are given in Table [2] 
The other T npq coefficients can be determined by applying the formula 

Tnpq — k n TQ p — n ,q-n- (213) 

In our classification, the shell model (1212l i is an N2-model. The corresponding function W(X, Y) is given by Eq. ( ID. Il l 
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Table 2: Numerical values of the Tq m coefficients )Frildll983l) . 



with real variables and coefficients 

W n (X,Y) = k n y i [ C\(X n - m -\Y n -\ — AX„-,„Y n+ i) + C2(X„-iY„- m -i — A m+l X n+m Y„ +m+ \) + Ci(X„+\Y n - m - A m X n+m+ iY n+m ), 

m>0 

(214) 

where C\, C2 and C3 depend on m only. Taking A = 2 and identifying the different terms in Eqs. (1212b and (1214b we 
find that the T npq coefficients must satisfy the relation 

Tn,n-m,n+\ 27 ''n,n—m—\,n—\ +2 7 \,n+m,n+m+\ ~ (215) 

which corresponds to energy conservation. The relation corresponding to enstrophy conservation is 

T n,n-m,n+\ 2 7 ' n s n—m—\,n—\ 2 T n n + m n + m +\ — 0. (216) 

Finally, from Eqs. ( 1213b . ( 1215b and (12161 l. we have 

Tn,n-m,n+\ — ^ri^O—m,! (217) 

T n ,n-m-\,n-\ = (1 _ 2 2m )(2 2 '" '-2) 1 k„To- m ,l (218) 

= 3 (2 -'"-2 2+m )- 1 A:„ro,- m , I . (219) 



Putting m = 1 in Eq. (1212b gives an L2-model with local interactions only. iFrikl (119831) found that such local inter- 
actions provide about 35 % of the total energy transfer and 25% of the total enstrophy transfer, showing the relative 
importance of non-local transfers. 

The generalization of model (1212b to MHD 2D turbulence, with conservation of total energy, cross helicity and 
square of magnetic potential (IFrikl 1 19841) is given by 

d,U n + vklU n = Y} T » 

,n-m-\,n- 

\(U n - m -iU„-\ - B 

n—m—l B„-i) + Tn,n-m,n+ 1 ( ^n-m U n + { B n - m B n+ \ ) 

;j;>0 

+ T \i,n+m,ri+m+l{U n+mU 'n+m+l ~ ^n+m^n+m+l)] + (220) 
,n-m-\,n- \(Un-m-\B„-\ ~ B 

n-w-l U„-i) + M„ \{U„- m B„ + \ - B n - m U n+ i) 
+ M nn+mn+m+ \(U n+m B n+m+ \ — B n+m U n+m+ i)], (221) 

with 

4 1 1 

^n,n— m— l,n— 1 — "J "^phn ^V* — m— l> n— 1' ^n,n—m,n+1 — ^ 2~^ m ^ n,n ~ m ' n+ ^' ^n,n+m,n+m+\ — '^Im^l jT^>yi+Hi,n-f-m+i(222) 

The solutions of Eqs. ( I220B22 It reproduce the expected properties of 2D MHD turbulence: a direct kinetic energy 
cascade and growth of enstrophy. In addition, the magnetic energy does not grow, meaning that a 2D dynamo is 
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impossible. However, there is an inverse cascade of the square of the potential vector. This implies that the magnetic 
energy spectrum becomes steeper in time, with the energy concentrated at the largest scales. The role of non-local 
interactions was not studied in this work. 

The hierarchical approach ab ove has been applied to various problems: passive scalar in 2D turbulence and 2D 
turbulent convection ( Frik , 1986 ). quasi- 2D convective turbulenc e in a thin vertical layer ( Barannikov et al. . 1 1988b . 
quasi-2D turbulent convection in a layer ( Aristov and FrikI 1990affb ) or in a rotating system ( Aristov and Frik , Il989 ). 

It is worth mentioning the study by Aristov and FrikI ( 1988b in which the hierarchical approach was used to model 
quasi-2D turbulent flow in a thin layer of an electrically conducting fluid heated from below. The layer was rotated, 
between two solid horizontal boundaries in a vertical applied magnetic field. For strong rotation and a weak magnetic 
field, the dimensionless quasi-2D equations for the 2D-fluctuations of the velocity u, magnetic field b and temperature 
9 are given by 



d,u + (u ■ V)u - (b ■ V)b + (t ■ V)t = - Vp + vV 2 u - yt/u, 

d,b + (u ■ V)b - (b ■ V)u = i]V 2 b, 

8,8 + (u ■ V)6» = kV 2 9, 

V u = 0, V b = 0, 



(223) 
(224) 
(225) 
(226) 



where t = (d x 8, 8 y 8} is the temperature gradient, fi characterizes the viscous friction at the horizontal boundaries, and 
v, if and k the viscosity, magnetic and thermal diffusivities. The set of ideal quadratic invariants is then 



(U 2 ) + <b 2 > - <T 2 >, 



<U b) - <U T>, 



h = <a 2 >, 



h = (e 2 ), 



(227) 



where a is the magnetic potential. In the limit 9 — > this set of invariants reduces to total energy, cross helicity and 



the square of magnetic potential as expected in 2D MHD turbulence. However, because in its general form Eq. ( 1227b 
includes two subtractions in I\ and h, a variety of scenarii are possible. For example, both (b 2 ) and (r 2 ) may grow 
simultaneously without changing I\ . 

The shell model corresponding to Eqs. ( 1223112251 is described by 



d t U„ = J] T npq [U p U q - B p B q + 9 p 9 q ] - vklU„ - fiU n , 

p.q 

d,B n = J] M npq [U p B q - B p Ug] - rjklB n , 

PA 

d,8„ = J] M npq [Up9 q - 9 p U q ] - Kk%. 



(228) 
(229) 
(230) 



p-q 



Some solutions are shown in Fig. [15] for different initial conditions. Only the bottom-right figure corresponds to a 
non-ideal case. The three others correspond to v = rj — k — 0. If at t = the temperature gradient is weak, top-left, the 
sum of kinetic and magnetic energies is constant, as it should be in ideal MHD. If at t — both the magnetic energy 
and temperature gradients are weak, top-right, they can grow without limit. If at t = the magnetic field is weak, 
bottom-left, the kinetic energy and temperature gradients grow simultaneously. On bottom-right, it is shown that even 
with non-zero dissipations, growth of the three quantities is still possible. In this case the 2D anti-dynamo theorem 
does not appl y. Because the temperature gradient cannot reach an infinite value, growth will eventually saturate. 
Originally |Zirmrj Q 198 lb introduced the hierarchical model for 3D HD turbulence, with 



u(t,x,y,z) = ^ U nNv (t)u nNv (r-r„ N ), 



(231) 



nNv 



where the third index v defines the orientation in space of the vortex, and can be equal to 1,2 or 3, denoting one of 
the three perpendicular directions of the 3D space. The shells are now defined in 3D Fourier space, and the base of 
function in the 3D space is 



u«jvv(r ~ r «/v) = C(s x e v ) (sin 2s - 2s cos 2s - sin s + s cos s) / s 



(232) 
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Figure 15: Time evolution of kinetic energy (thick-red), magnetic energy (thin-green) and square of temperature gradients (dashed-blue) obtained 
from the shell model Eqs 122812301 , corresponding t o quasi-2D turbulence in a thin rotating layer of conducting fluids, for four different initial 
conditions and ideal or non-ideal cases. Adapted from lAristov and Frik 1 1988). 



The density of vortices now increases with n as 2 



3« 



From this base of function a new shell model for 3D turbulence can be constructed following the approach of Frik 



(119831) b ut with helicity as second ideal invariant. To our knowledge this remains to be done. An attempt has been 



made by Shaidurova (1987) but enstrophy was still use d as second ideal invariant 



Also using the base of function given by Eq. (1232t . lZimin and Hussainl (11995b obtained an Nl -model given by 



d,U„ + vk\u n = k n £ A m [A- 5m/2 U n - m U„ - A- 3m/2 U 2 n+m ]. (233) 



m>~\ 

where A_i = 0.387, A m >j = 2.19 and Ao = 0. Only kinetic energy is conserved in this model. However, one 
remarkable feature is that the infra-red scaling laws are well reproduced. A generalization of this model to MHD is 
proposed in Sec. 13.41 

3.3.2. Non-local version of the Sabra model 



The non-local version of the Sabra model for 3D MHD turbulence has been introduced bv lPlunian and Stepanov 



(120071) . It corresponds to 



W„(X, Y) - — ^ A m [ A"'(l + A)(X* +m Y„ + ,„ + i + X„ +m+ \Y* +m ) + (-l) m+1 (X^ +m F n+m+ i - X„ +m+ iY* +m ) 
-A(l - (-A) '" l )(X*_ m Y n+1 + X„ +1 Y*_ m ) + (X*_ m Y„ + i -X„ +1 Y*_ m ) 

+ x n 

—m—\ *n—\ ) ] , (234) 

leading to the system 

d t U n = i ^ A„, [ (k„ +m + k„ +m+ i)(U* +m U„ +m+ i - B* n+m B n+m+ \) - (k n +\ + (-l) m k„- m )(U*_ m U„ + i - B* n _ m B n+ \) 

+ (k n -i + (-l) m+ Vm-i)(t/ n -i C/»-m-i - B„-lB„- m -l) ] - vk\U n + F„, (235) 

^A m [(-iy" +i (f/: +m B„ +m+ i - B* n+m U n+m+ i) 

m>\ 

+ U*_ m B n+ \ - B* n _ m U„ + i + U n -\B n - m -\ - B n -\U n - m -\ ] - T]klB„, (236) 
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where A m is now some arbitrary parameter depending on m. This model conserves total energy, cross helicity and 
magnetic helicity, the latter being defined by Eq. ( 11951 ). The choice of interacting triads is not arbitrary and corresponds 
to all possible triads (Sec. 13.1. 21 . 

The Sabra model for 3D MHD turbulence corresponds to m - 1 and Ai = (A + A 2 )~ l , 



d,U„ - ik„ [ U* +l U„ + 2 - B* n+l B n+ 2 -p— (U*_ l U n+ \ - B* 1 _ l B„ +i ) 

+ -j ( U n . i U n - 2 - B„- iB„_ 2 ) ] - vkl U„ + F„ , 



d t B„ = 



A + A 2 



[ - B* n+l U n+ 2 + U*_ l B n+ i - B*_jC/„ + i + U n -\B n -2 - B n -\U n -2 ] ~ ^k 2 B n . 



(237) 



(238) 



For B — 0, the HD model introduced bv lL'vov et al.l d 1998b is refound. It conserves kinetic energy and helicity, 
the latter being defined by Eq. ( 1187b . 

In contrast with its local version, the sy stem of Eqs. d235H236b is not uniquely defined as it depends on the 



Plunian and Stepanovl (12007b suggested that A m = A y(m - 1] /A(A + 1) with y < 0. The Sabra model 



parameter A 

corresponds to y — > —oo 

Assuming isotropy, iPlunian and Stepanovl ( 12007b estimated y that accounts for the number of all possible triads 
between three shells. They found y - -7/2 for N2-models and y = -5/2 for Nl-models, the latter being consistent 



Zimin and Hussainl ( 1995b 



with the helical model derived by 

For free decaying turbulence, Plunian and Stepanovl ( 2007 ) found that y does not significantly change the slopes 
of the kinetic and magnetic energy spectra in the inertial range, both compare well with Kolmogorov scaling \U„\ 2 oc 
\B n \ 2 oc k^ 2 ^ 3 (corresponding to a k~ 5 ^ 3 energy density spectrum). On the other hand, y has a strong effect on the slopes 
of the infra-red spectrum, as shown in Fig. [16] It corresponds to the following dependency 



Block: 2 ?* 2 . 



(239) 



Different values of y may account for the presence of different infrare d mechanisms in real 3D HP a nd MHD 
turb ulence. We note that y = -5/2 is in agreement with the arguments of iBatchelor and Townsend 
and Pouauet et al.l d 1976b fo r MHD. However, other values of y could also be satisfactory dSaffman . 



1967 



1948) for HD, 



Lesieur 



1997FFournieret all 1982 ). We note that the Sabra model given by Eqs. (1237H238t . cannot give a realistic slope 



for the infra-red spectrum, emphasizing the importance of including non-local interactions. Finally, for forced MHD 
turbulence (dynamo action) there are strong phenomenological arguments for taking y — -1 as explained in Sec. 14.1. j] 



3.4. Helical models 

The kinetic and magnetic helicities given by Eqs. ( 1187b and (1195b . have the particularity of being stro ngly corre- 
lated to, respectively, the kinetic and magnetic energies. In particular, as noted bv lBiferale and Kern d 1995b for the HD 
case, the kinetic helicity "presents an asymmetry between odd and even shells that does not have any counterpart in 
physic al flows". The same remark applies to magnetic helicity. In order to circumvent this problem. lBiferale and Ken 
(11995 ) introduced another type of shell mo del based on the decompos ition into helical modes in the manner described 
bv lWaleffel d 1992b . At about the same time. lZimin and Hussainl d 1995b also introduced a helical, non-local shell model. 
The generalization of these helical shell models to MHD is given below. 

In order to account for the helical decomposition of the velocity and the magnetic field given in Eqs. (1861187b . we 
double the number of variables per shell, U~,B+ and B~. The kinetic and magnetic energies, the cross helicity, 
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-4 -2 2 4 6 



-4 -2 2 4 6 



Figure 16: Free decaying MHD turbulence for v = 10~ 6 and Pm = 1, with (top-left) local model, (top-right) y = -2.5, (bottom-left) y = —1.5 and 
(bottom-right) y = — 1. Grey corresponds t o kinetic energy, red to magnetic energy. Increasing time is indicated with spots of lighter intensities. 
Adapted from Plunian and Stepanov] i2007t) . 



H° 



and the kinetic and magnetic helicities are now 

e u = ^K |2 + |f/ « |2 ) (240) 

eB = ^ZK |2 + |B « |2 ) (241) 

n 

HC = \Yl( U « B « ' + U n B n "+CC.) (242) 

/J 

^^(|f/,tl 2 -|f/,7l 2 ) (243) 

n 

2^'Kl 2 -|B,Tl 2 )- (244) 



The enstrophy and square potential become 

S = ^fc 2 (|t/,!l 2 + |f/,7l 2 ) (245) 

A = ^2^ 2 K |2 + |B » |2 )- (246) 

The general formalism derived in Sec. 13.1.21 and given by Eqs. ( II 45111 46b still applies, defining any vector X as 

X = (X+, X\, ■ ■ ■ , X+, X-, X 2 , ■ ■ ■ , X N ). (247) 
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Applying this definition to U, B and W, we can rewrite Eqs. ( 1145111461 in the form 



and Eq. ( TBTl as 



d t U$ = W,f(U,U)-W,f(B,B)-v^f/,t +K> 
d t B± = W^(U J B)-W^(B,U)-» 7 ^B* 

J] W-„ + (X,Y)C + W-{X,Y)Y-* + c.c. = 0. 



(248) 
(249) 

(250) 



In principle shell-to-shell and mode-to-mode energy transfer and flux can also be obtained using the definition 
of X given in Eq. (1247b . Transfer involves not only the velocity and the magnetic field but also both helical modes. 
Instead of just four, they are now sixteen possible transfers, t/ ± -to-t/ ± , U ± -to-B ± , B ± -to-U ± and B ± -to-B ± . 

We note that in each shell «, only the helical modes t/* and Z?* are defined, from which the energies and helicities 
are calculated. However if U* and B^ are taken to be real, then new complex variables U„ and B n can be defined 



E u 



with energies and helicities 

n 

= ^I>«i 2 

n 

= I J] (U n B n * +c.c.) 

n 
n 
n 

n 
n 

The equivalence between Eqs. (1240H244t and Eqs. d245H246t is straightforward. 



enstrophy and square potential 



(251) 

(252) 
(253) 
(254) 
(255) 
(256) 

(257) 
(258) 



We call HI -model a first-neighbor model for which magnetic helicity as defined in Eq. ( 1244b or ( 1256b is conserved. 
Similarily we call H2-model a two-first neighbor model with the same quadratic invariant. In Sec. 13.4.11 we present 
the general expression of Hl-models derived in terms of complex variables U„ and B n (real {/* and B„). In Sec. 13. 4. 21 
several H2-models in terms of complex variables U* and B„ are presented. 

3.4.1. Hl-models (helical, two feet in the same shell) 



A he lical version of th e HD m odel elaborated by| Zimin and Hussain ( 1995) was first used to study 3D HD turbu- 
lence bv lMelanderl (119971) (see also M elander and Fabiionasl (120021 2003, 2006)). The model is based on real variables 
U* and U~ and is both helical and non-local. Equations were given for variables S„ = (U* + U~) and D„ = (U* - U~) 



43 



and after appropriate renormalization we obtain 



d,S „ = k n Yj a »> {S nS n-m - r m D n D n . m - A m S ; +m + A m D 2 n+m ) - vk;S „ + F s n . 

m>-l 
m#0 

d,D„ = k„ A». (F m S n D n . m - D n S n - m ) - vk 2 n D n + F°. 



(259) 



m>-\ 
m*0 



where A m is again an arbitrary parameter depending on m as in model (1234b . Using a hierarchical approach, Zimin 
found that A m = A- 5[m+1/2]/2 . However, any other function of m does not change the c onservation of kin etic energy and 
helicity. This model was solved mainly for free decaying turbulence (F^ = F® = 0). iMelanderl ( 1 19971) found that any 
solution starting from a point close to the Kolmogorov stable manifold in phase space, misses equilibrium and enters 
a chaotic regime characterized by an exponential growth in helicity as the solution diverges from the equilibrium. 
Erratic flu ctuations of helicity and their d issipation imply intermittent energy decay. Early enstrophy divergence wa s 
studied bv lMelander and Fabiionasld2002h and the transient laws were investigated bv lMelander and Fabiionasld2003l) . 
With 

U n =iS n ±D n , (260) 
Stepanov et al. I d2009h rewrote the local version of model ( 1259b in a complex form. After appropriate renormalisation 



d,U„ + vklU n = ik n [Ul_ x + U* n _ x 2 + AU* n (U n+1 - U* n+l ) - A 2 U„(U n+l + U* +1 )] 

- cik n [u„(u n ^ + u* n _,) + A.u* n {u* n _, - c/ n _!) - a\u 2 1+x + u; + 2 )] + F n , 



(261) 



where c is a free param eter. The kinetic energ y and helicity defined by Eqs. (1252b and (1255b are ideally conserved 
whatever the va lue of c. Stepanov et al. ( 20091) chose c = A~ 5 / 2 in order that model (1261b becomes equivalent to the 
local part of the lZimin and Hussainl (11995b model. In terms of S„ and D„, the kinetic energy and helicity are given by 



11 

H u = ±2^k„S„D„. 
The version of model ( 12611 ) corresponding to 3D MHD turbulence is given by dMizeva et al. 1 12009I) 



(262) 
(263) 



W„(X,Y) 



ik„ (X„_iF„_i + X^Y^^) - AX*Y* +l — ^-(^1^,1+1 + X„ + iY„) 

~~^P^n-\Yn-\ — ^n-\Y„_{) + AX n Y n+ \ ^"(^n^n+1 +^n+l^n) 



ick„ 



(Xn-iY n +X„Y„-\) + AX* 1 Y*_ l - A 2 (X n+ \Y n+ \ +X* +l Y* +1 ) 



+ 2~(XnYn-i +^ r „-i^ / «) ^X n Y„-i + —(X n+l Y„ + i — X n+ iY l1+l ) 



(264) 



for which the magnetic helicity ( 1256b is conserved in addition to total energy and c ross helicity. This model was used 
to study MHD helical t urbulence dMizeva et all 2009t iFrick and Stepanov! l2Q10t) . including global rotation and an 
external magnetic field ( Plunian and Stepanovl 2010 ). The results will be detailed in Sec. 14.4. ll 

3.4.2. H2-models (helical, two feet in two neighboring shells, third foot in a smaller shell) 



Expanding the velocity Fourier modes onto a base of polarised helical waves as described in Sec. 12.2.51 Biferale and Kerr 



d 19951) introduced four independent helical shell models for 3D HD turbulence. Following iLessinnes et al.l d2011l) 
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these can be written in a compact form 



d t U$ + vk 2 n Ut = ik n [( Sl A - s 2 X 2 )U*l\U*% + (s 2 A - A- l )U^\U^- + (A~ 2 - si^U^U??]' + F± 



(265) 



where {/* are complex variables. The four independent models correspond to (s\, s 2 ) £ {(1, 1), (— 1, 1), (1, -1), (-1, -1)}. 
The kinetic energy and helicity defined by Eqs. ( 1240b and (12431 ) are ideally conserved. 



For (si, s 2 ) = (-1, 1) Eq. ( 12651 ) consists of two independent sets of equations forU = (t/f, U- , Ut , ■ ■ ■ , U 2 _ v U 2 ) 
and U = (t/ r, Uif, U^, ■ ■ ■ , UZ. 1, UZ S), each of which corresponds to the GOY model given by Eq. ( 1185b with 
s = 1 - A' 1 dBiferale and KentTl995b . 



• For (ii, 52) = (1, 1) Eq. (1265b consists of two independent sets of equations for U = (U^, U 2 , • ■ ■ , Ui_ v U 2r ) 
and U - (U7, U 2 , \ ' ■ 1 ^ 2n _ v ^2 ) eacn °^ w hich conserving separately a positive-definite quantity. This lead 



Benzi et alj dl996al) to conclude that such a model "is equivalent to two uncorrected GOY models for 2D 



turbulence". 



For (s\, so) = (-1,-1) after lBenzi et al.l d 1996al) the model "may present a significant backward energy transfer, 
which leads to possibly strong deviations from the Kolmogorov scaling". 

For (^1,^2 ) = (1, -1) the mo del shows intermittent statistics in excellent agreement with the Navier-Stokes 
equations (IBenzi et al.Ul996al) . It is given by 



d t Uf, + vklut = i [k n+l (l + A)U± +l U* +2 - k„(A + A~ l )U^ + ft^A" 1 - Wn-2 U n-i\* + F t- ( 266 ) 



Another version of this model has also been formulated bv lChen et al.l d2003l) . in the spirit of the Sabra model 
d t U$ + vklUt = i [k„ +l (l + X){U$ +1 )*Ul+ 2 - k n {A + A-'W^TU^ - k,^(A- ] - DU^U^] + Ff x . (267) 

The interest of helical rather than non-helical models has been clearly demons trated by the study of the helicity flux 
in 3D HD turbulence. Using the GOY model, Ditlevsen and Giulianil ( 2001allb ) found a range of scales, within the 
inertia l range, for which the helicity spectral flux is exponentially divergent in k due to dissipation (see also lDitlevsen 
(1201 lb ). This implies that even if the flow is helical at large scales, there is a subdomain of the inertial range for which 
the flow is non-helical. This picture contrasts strongly with that arising due to the energy cascade occurring all along 
the inertial range for which dissipation takes place only at scales smaller than the Kolmogorov scale. The physical 
reason why helicity and e nergy do not dissipate on the sa me scale was not clear until helical models came into being. 

First, using both a HI (IStepanov et al. ■ l2009h and H2 dChen et al.Ll2003l) models, the helicity spectral flux has been 
found to be constant all along the inertial range without the divergence in k previously found with the GOY model. 
Apparently both types of model, non-helical and helical, were giving contradictory res ults. 

Then using an H2-model corresponding to (s\ , S2) = (1,-1 uLessinnes et al.l (12.01 II) calculated the helicity spectral 
flux for both + and - helical modes and found a subdomain of the inertial range in which both helicity flux spectra 
diverge in k due to dissipation in agreement with the GOY model. However, as both fluxes have opposite signs they 
also found that their sum, which is the flux of total helicity, has a constant spectrum all along the inertial range. In the 
GOY model two helical modes cannot be present in the same shell (see the case (si, s 2 ) = (— 1, 1) above) and therefore 
they cannot be summed. This clearly shows the superiority of helical shell models in dealing with kinetic helicity in 
3D HD turbulence, and presumably with magnetic helicity in 3D MHD turbulence. 

The helical shell model of MHD turbulence elaborated bv lLessinnes et al. I d2009bh is a linear combination of the 
four previous submodels (si, s 2 ) - (±1, +1), and includes an estimate of the weighting to be given to each submodel. 
In fact it is not clear whether such a combination is appropriate, as only the model for (si,S2) = (1,-1) gives a relevant 
helicity spectral flux for 3D HD turbulence. Instead it might be more appropriate to consider the MHD version of 
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models ( 12661 ) or ( 1267b . Their respective expressions in terms of W ± (X, Y) are 

wj(x, y) = i (fe n+1 (i + i)(^ +1 c 2 + Ci^ +2 ) - - w 

and 

w,f(x, Y) = i [k n+1 (i + A)(x-yj +2 + O^J + 2) - fVdC 2 - Ci*^) 

- wXUK-i + - 4-iC 2 C - C 2 ^i)) ■ ( 269 ) 

In both models total energy, cross helicity and magnetic helicity are ideally conserved. 



4. Applications of MHD shell models 

With the exception of models ( 1268b and ( 1269b which are introduced here for the first time, the other shell models 
introduced in Sec.[3]have been used to investigate the statistical properties of forced and free-decaying MHD turbu- 
lence (Sec. 14. H and 14.2b . They have also been the starting point for more complex models, like multi-scale dynamo 
models, Alfen-wave models or Hall-effect models (Sec. 14.31 I4.4l and 14.5b . 



4.1. Forced MHD turbulence 
4.1.1. Forcing 

In forced MHD turbulence, choosing an appropriate forcing F is of course essential. It is generally applied at the 
largest scales of the system. The two main features of a forcing are its time dependency and the range of scales over 
which it is applied. 

With a stationary complex forcing F applied on one scale, the injection rate of two real quantities, e.g. energy e 
and cross helicity x, can be controled. However, with such a stationary forcing the solutions of Eqs. ([12 1H 122b may 
reach an unphysical sup ercorrelated state U„ = +B„, with no fluxes and no inertial ranges (IGiuliani and Carbone . 



19981 : iFrick et al.L 120001) . Instead using a forcing F„ = f n e l<f,M , where <p n (f) is a random phase and f„ is real, leads 
to long lasting dynamics of U and B. However, using such random phase, the injection rate of only one quantity 
instead of two is then controled. For example, taking <p„(t) constant during time interval t c , but changing random! 



from one time interval to the next, the injection rate of energy becomes e * f„t c (IPlunian and Stepanovl 120071 1201 



my 
10) 



and the cross helicity injection is not controled anymore. I n order to control the injection rate of two quantities again, 



the forcing must be split between two neighboring scales (IMizeva et al 
following set of equations 



2009). In that case the forcing satisfies the 



j bf+I 

- ^ u;f„ + u„f; = e 

n=np 



(270) 
(271) 



Mizevaet al. (2009) used their helical shell model given in Eq. (1264b with Pm = 1 to study the influence of cross 



helicity on MHD turbulence. They found that the total energy of the system increases with the ratio x/e, following 
the scaling law E oc (1 -^/e) -4 ^ 3 , as indicated by the full curve in Fig.[17](left panel). This corresponds to an increase 
in the non-linear transfer time by a factory. In other words, the energy transfer in the inertial range is depleted, 
implying energy accumulation in the large scales, and a steeper spectrum as shown in Fig.[T7](right panel). In the limit 
X/e — > 1 the velocity and magnetic fields are correlated, the non-linear terms in the MHD equations are canceled, and 
the energy cascade transfer is blocked. 
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Figure 17: Left panel: Total energy of stationary forced MHD turbulence versus the relative cross helicity injection rate x/e. The time evolution of 
energy is shown in the inset, for^ = (thick line), x — 0.3 (grey line), and^ = 0.6(thin line). Right panel: Total energ y spectrum norm alized by 
iT 2 ' 3 for three different injection rates of cross helicity. The straight line corresponds to E(k) oc it -1 - 9 . Adapted from lMizeva et al.l j20O9t) . 



One may also want to c ontrol the injection rate of kin etic helicity. Depending on the model used t his might require 



an add itional forcing shell (IStepanov and Plunianl 120061) or forcing on two different helical modes dLessinnes et al 



2009b) 



One could also envisage applying some forcing G in the time evolution equation for B. Such magnetic forcing is 
naturally implemented in multi-scale dynamo models to accoun t for the back-reaction of the large-scale fields onto 
small-scale turbulence ( Frick et al. , 20061 : Nigro and Veltri , 201 1 ), and will be detailed in Sec. 14.31 Note that magnetic 
forcing can occur without injecting magnetic energy. Let us take the example of constant injection rates for kinetic 
energy e, cross helicity x ar, d magnetic-helicity in the absence of magnetic energy injection. Consequently we 
can take F and G to be constant in time (no random phase) within shell rif . With the magnetic helicity defined by 



Eq. ( 12561 ). F np and G„ F must satisfy 



Uu* nF F nF + U nF F* nF ) = 6, 



K f G "f + B n F G*„ F 



l - (b: f f„ f + b„ f f: f + u; f g„ f + u„ f g: f ) = x , 

ik n] {K F G*„ F - B„ F G nF J = 



(272) 
(273) 
(274) 
(275) 



The results obtained using such forcing applied at k„ F = 1 are given in Fig. [18] An inverse magnetic helicity cascade 
is found corresponding to a negative flux of magnetic helicity for k < k„ F (right panel). 



4.1.2. Small-scale dynamo action 

We speak of dynamo action when the magnetic energy is sustained by the motion of an electroconducting fluid. 
For the dynamo to work the magnetic Reynolds number, e.g. at the forcing scale of the system, must be higher 
than some threshold value Rm c . Experiments and simulations show that Rtry depends on the presence of a mean 



flow and its characteristics, such as geometry and time dependency (IPevrot et al 



201 II). the electromagnetic boundaries, such as conductivity a nd permeability dFrick et al 



2007 



2008; Ponty and Plunian, 



2002; Dobler etal., 2003 



Avalos-Zuniga et al.ll2003tlAvalos-Zuniga and Plunianl 120051) . the presence of an external magnetic field, and global 



rotation. 

The dynamo threshold Rm f can also be calculated assuming that dynamo action is produced only by the turbu- 
lence. Then dynamo is reached only if a sufficient amount of power e c is injected into the fluid, corresponding to 
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Figure 18: Energy, cross helicity, magnetic helicity spectra (left panel) and fluxes (right panel) of forced MHD turbulence with constant injected 
helicities (e = \,x = 0-3, £ = 0.4). 



Rm c = ej^tp jr\. In Fig.[19]fhe dynamo threshold Rm e is plotted versus Re ( Buchlin. E. . 201 1 ). The curves were 
obtained using the N2-model given by Eqs. (1235H236t and three degrees of non-locality: local y = -co (solid line), 
weakly non-local y = -5/2 (dashed line) and strongly non-local y — -1 (dotted l i ne), w here the parameter y is 
defined in Sec. 13.3.21 The results are similar to the ones obtained by IStepanov et al.l (12006) using a 3D MHD GOY 
model, suggesting that the results are sensitiv e neither to the mod el (GOY or Sabra) nor to non-locality. The dia- 
monds correspond to Rm f obtained from DNS (Iska kov et al.l 120071). The overall picture is that the dynamo threshold 
first increases with Re, as found from DN S by Schekochihin et alT ( 2004 ). and then decreases to reach a plateau for 
Re > 10 4 . Similar results were obtained bv lPontv et al. I d2005h using Large Eddy Simulation technics 
Usin g a helical shell model derived from a combination of the four H2-models given by Eq 



Lessinnes et al 



(12009bl) compared the effect of maximum helical and non-helical forcing. They found that the dynamo threshold is 
lowered by at least a factor 1.5 when using a maximum helical forcing. 
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Figure 19: Comparison of the dynamo threshold Rm c versus Re, obtai ned from the N2-mo del given by Eqs. <235I2361 with different de eree of 
non-locality (the three curves) and from a DNS shown by the diamonds (Iskakov et alll2007l) . From Buchlin ~Ej 1201 lh . 



Starting with a statistically stationary turbulent flow, the time evolution of the magnetic energy for Rm > Rm r 
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can be split in two regimes, the kinematic regime during which the magnetic energy grows exponentially, and the 
saturated regime starting when the magnetic energy reaches some statistical stationary level. Other dynamics like on- 
off intermittency or chaotic reversals, which involve large-scale fields, are not captured with shell models unless large- 
scale equations are added (see Sec. 14.31 ). In Fig. [20] two examples are shown for Pm «: 1 (left-panel) and Pm » 1 
(right-panel), and Rm » Rm f in both cases. The magnetic energy spectrum (in red) grows with time. The end of the 
kinematic regime corresponds to the point at which the flow becomes sensitive to the magnetic field via the Lorentz 
forces. For Pm <K 1 (left-panel), a depletion of the kinetic energy spectrum occurs due to the equipartition between 
magnetic and kinetic energies. A significant part of the injected power e is then lost through magnetic dissipation, 
explaining why the viscous scale l v becomes larger (k v becomes smaller). For Pm » 1 (right-panel) refilling of the 
kinetic energy spectrum occurs in a range of scales between l n and l v , due to direct transfer from magnetic to kinetic 
energy. The results of Fig.[20]have been obtained with, again, the N2-model given by Eqs. ( 123511236b . with y — -5/2 
for Pm = 1(T 3 , and y = -1 for Pm = 10 4 . 




Figure 20: Kinetic (black) and magnetic (red) energy spectra at different times for Pm = 10 3 (left-panel) and Pm = 10 4 (right-panel) . The time 
interval between two samples is constant. Along time, the curves evolution is given by the arrows. Adapted from Stepanov and Plunian (2008). 

For Pm <k 1 the results are consistent with dynamo action occurring on a scale fci of the same order of magnitude 
as the magnetic dissipation scale k n l , and a kinematic growthrate T^n such that 



^kin 



cc ^/y 3 / 4 . 



(276) 



This corresponds to the dynamo action being produced by local energy transfers, the main physical mechanism being 
the action of fl ow shear against magnetic d issipation. This is also consistent with fast and small-scale kinematic 



dynam o action ( Childre ss and Gilbert! 1 19951) . Taking other values of y e [-co, -0.5] does not significantly change the 



results ( Plunian and Stepanov , 2007h . 

For Pm » 1 the parameter y is taken to be equal to -1 in order to have dynamo action occurring at the viscous scale 



^kin ~ k v , where the flow shear is maximum dSchekochihin et al 
regime is then characterized by 



2002; Schekochihin et al., 2004). The kinematic 



r kin = k v u(k v ) cc e 1/2 v" 1/2 , 



cc e 1/4 v- 3/4 . 



(277) 



The model is non-local and the kinematic dynamo is again fast and small-scale. 

The transient regime leading towards saturation has also been studied in detail for Pm <k 1 and Pm » 1 . An 
inverse cascade of magnetic energy occurs towards large scales. In the saturated state the ratio of magnetic to kinetic 
energy is found to be larger than unity (about 1.5). Such super-equipartition is spread over the whole inertial range. 

In such a saturated state, characterized by a statistical stationarity for both kinetic and magnetic energies, the 
time-average of the growthrate of the magnetic energy is thus equal to zero. Does it imply that the saturated flow 
is u nable to produce dynamo ac tion anymore and can just compensating for dissipation? To investigate this ques- 



tion 



Cattaneo and Tobias d2009h introduced a passive vector c, satisfying a third equation identical to the induction 
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equation: 



{d,-rjV 2 )c = -(u ■ V)c + (c ■ V)u, V ■ c = 0. 



(278) 



Here the vector c is passive in the sense that it does not back react onto the flow. In other words, the system of 
Eqs. ([T]|2]i remains unchanged. Cattaneo and Tobias d2009h solved the problem for Pm = 1 using two approaches, a 
DNS and the MHD GOY model given by Eqs. j 193H 194-b . Both calculations lead to the same conclusion summarized 
in Fig. [2J] for the GOY results. After a saturated regime has been reached (top figure) the additional Eq. ( 12781 l is 
switched on. The passive vector energy is found to grow exponentially (bottom), showing that the saturated flow 
keeps its small-scale dynamo characteristics, presumably due to its chaotic dynamics. 




Figure 21: Top panel: Time series for the kinetic energy (solid line) and magnetic energy (dashed line) showing the dynamo evolving to a saturated 
state (v = rj = 10~ 6 ). Bottom panel: Time series for the kinetic energy (solid line), magnetic energy (dashed line) and the energy of the passive 
field (dot-dashed line). Adapted from Cattaneo and Tobias (2009). 



Using again the MHD GOY model lLessinnes et al.l (I2009al) calculated a few of the fluxes introduced in Sec. 13.1.31 
for the saturated regime. These are shown in Fig.[22]for v = 10~ 9 and Pm = 1CT 3 , together with the kinetic, magnetic 
and total energy spectra. Rewriting Eqs. (1168H1691 for « > , in the form 



(279) 
(280) 



we define 



11 all 




+ n£, 


(281) 




ll u < + ll u> 




(282) 




_ 11 all + 11 all 




(283) 




ll B< + LL B> 




(284) 



We stress again that fluxes are always «-dependent though not explicitely appearing in the above notation (see 
Sec. 13.1.3b . The blue dot-dashed curve representing in the right panel of Fig. [22] is constant in the range 
kf — 1 and log 10 A: v « 6 where n^,' = e as 0.98. Such a constant value necessarily leads 



k 6 [k F ,k v ] with log I0 
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to a Kolmogorov spectral slope, as shown in the left panel of Fig. l22lFI. For scales smaller than the viscous scale 
(k > k v ), II^j drops by a value equal to the viscous dissipation e y ~ 0.34, leading to fl^ = e - e v — e n ~ 0.64. Note 
that the curves given by the blue-dashed and blue-dotted lines are symmetric about a horizontal line for k e [kf , k v ], 
implying that the relation = + is indeed satisfied in this range of scales. The blue-solid curve is above 
the blue-dot-dashed curve for k e [kp, k^], with log I0 k n « 3.5. This indicates that FPl > for this range of scales. 
Finally, the growth in the green-solid curve for k e [kf, k v ] with Yl^ > shows that dynamo action occurs at all scales 
in this range. 




As depicted in Fig.[22]the power e injected into the flow at the forcing scale is dissipated at smaller scales by Joule 
and viscous dissipation, according to 

€ = €„ + £y. (285) 



On the other hand, there is no theoretical argument predicting how e is distributed between e n and e Y . IPlunian and Stepanov 
(2010) studied this problem for different values of Pm, using the helical shell model given by Eq. (12641 1. In Fig. [23] 
the ratio r = e n /e v versus Pm is plotted for different values of v and rj. For a given value of v, in the limit Pm — > 0, 
r] +oo and dynamo action becomes impossibe implying r — » 0. For Pm = 1 both kinetic and magnetic spectra are 
identical, implying e y = e n = e/2, and consequently r = 1. For a given v, there is always an intermediate value of Pm 
for which r reaches a maximum. This is related to a super-equipartition state in which the magnetic energy is higher 
than the kinetic energy at large scales. For typical values of MHD turbulence, v < 10~ 8 and Pm as 10~ 5 , r > 10. 



4.2. Free-decaying MHD turbulence 

The decay properties of MHD turbulence have important applications in astrophysics. For example, the compar- 
ison of lifetime of a magnetic astrophysical object with the decay time of its magnetic energy can help us decide 
whether the magnetic field is of primordial origin, or produced e.g. by dynamo action. In problems of this nature 
the magnetic helicity is known to play a crucial role. The energy decay l aw satisfies E(t) oc r 1 or E(t) oc r 1 ^ 2 , 
depending on whether the magn etic helicity is respectively zero or maximal ( Biskamp and Mulleri 1999t Campanelr3, 



2004; Christensson et al. 



and Antonov et al 



2005). By solving the MHD GOY model given by Eqs. 



Antonov et all (2001) 



(1200 II) found that cross helicity a lso plays a crucial ro le. 

Using the helical shell model given by Eq. ( 12641 1. Frick and Stepanov ( 2010h explored the combined roles of cross 
helicity and magnetic helicity. Taking v = rj = 10~ 5 , they considered a set of 128 different initial conditions, with 
both kinetic and magnetic energies concentrated in the first shell n = 0, corresponding to ko = 1 (E¥ = E% « 1 
and t/„>o = B n> o = 0). A very small amount of magnetic and cross helicities QHq | < 10~ 4 and \H^\ < 10~ 4 ) was 
injected. In Fig. [24](left panel) the total energy versus time is plotted for the 128 initial conditions. At an early stage 



7 A -2/3 energy spectral slope corresponding to a -5/3 energy density spectral slope. 
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Pm 

Figure 23: Dissipation ratio versus Pm. The solid lines from right to left c orrespond to v = 1(T 5 1(T 6 1(T 7 , 1(T 8 . The dashed curves from bottom 
to top correspond to rj = 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 1/256. From lPlunian and Stepanovl feOlOh . 



(t < 30) we see that E(t) oc f°, with a e [-1,-1/2] (dashed lines), in good agreement with the DNS results and 
theoretical predictions. At a later time (? > 100) many of the curves decrease with an exponential decay, in accord 
with pure dissipation. On the right panel of Fig. [24] the ratio H c /E versus time is plotted. The curves on the left 
with an exponential decay correspond to the curves on the right with H c /E — > ±1. An H c /E as ±1 corresponds to 
U„ ~ +B„, implying again (see Sec. I4.1.U a depletion of energy transfer in the inertial range and the accumulation 
of energy in the largest scales where the exponential decay occurs. Note that such a final state is reached at rather 



long ti mes (t > 10 2 to 10 3 ), which may explain why this is not observed in DNS. As shown in iFrick and Stepanov 
l l2010l) . cross helicity is generated at the dissipation scale, and then slowly cascades backwards towards the large 
scales. This explains how a maximal cross helicity final state can be generated from a minimal cross helicity initial 
state (\HC/E \ < 10~ 4 ). Finally, in the left panel of Fig. \2A\ a few curves correspond to E(t) oc t for t > 100, with a 
limit H c IE + +1. The latter was observed previously bv Brandenburg etaL ( 1996) and is presumably relat ed to the 
concentration of magnetic helicity, rather than cross helicity, in the largest scales ( Frick and Stepanov , 2010h . 




Figure 24: Free decay of tot al energy E (left panel) and normalized cross helicity H c /E (right panel) for Pm = 1. The 128 realizations differ by 
their initial conditions. From Frick and Step anovl feOlOh . 



To complete the study of free-decaying MHD turbulence we have calculated, using the same helical shell model 
given by Eq. ( 1264b . the kinetic energy, magnetic energy, magnetic helicity and cross helicity spectra obtained at time 
t = 1000, with four different combinations of initial conditions: maximal or zero for both cross and magnetic helicity. 
The results are presented in Fig.l25lwith 
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(a) (H c /E)t=a - (H /E), = q = 0, the energy decays and helicities change their sign with k. 

(b) {H c /E) t= Q = and (H B /E) t =o = 1> the magnetic energy cascade is accompanied by a simultaneous inverse 
cascade of magnetic helicity towards the largest scales. This leads to the accumulation of energy at the largest 
scale (the energy spectra are extended to the left). 

(c) (H c /£),=() = 1 and (H B /E), = o = 0, the cross helicity blocks the energy cascade, leading to steeper spectra. 

(d) (H c /E), = o = (H b /E), = q = 1, both effects are combined with essentially weak non-linear energy transfers. 




Figure 25: Spectra of kinetic energy (blue squares), magnetic energy (magenta circles), cross helicity (brown triangles), magnetic helicity (green 
diamonds) obtained at time t = 1000, from free-decaying MHD turbulence (v = T] = 10~ 5 ) and normalized by <T 2 ' 3 . Each panel corresponds to 
different initial conditions with (a) H^_ Q = Hf =Q = 0, (b) = 0,Hf =Q = 1, (c) fl£ = 1,#* = (d) H^_ Q = Hf =Q = 1. Filled (empty) symbols 
correspond to positive (negative) values of helicities. 

4.3. Multi-scale dynamos 

Many MHD applications are primarily concerned with the properties of large-scale quantities, e.g. the chaotic 
reversals of the Earth's dipolar magnetic field, the dynamics of sunspots presumably due to underlying large-scale 
magnetic structures in the convection zone of the Sun and the zonal winds in Jupiter. Because of obvious numerical 
limitations it is not possible to solve the equations over all scales, rangering from the smallest dissipation scale to the 
size of the object itself. This is why we need a model for small scales (e.g. a shell model), so that only the large- 
scale equations remain to be solved (e.g. by DNS). Such splitting, between large and small scales, makes it crucial 
to understand to what extent they are linked, namely how do small-scale effects influence large-scale structures and 
vice-versa. Beyond numerical issues, this problem is at the heart of the physical understanding of HD and MHD 
turbulence. 
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4.3.1. Reynolds equations for MHD 

Both velocity and magnetic fields are split into large-scale and small-scale components 



u = u + u 



b = b + b', 



(286) 



e.g. after applying a large-scale filter. On similarly splitting Eqs. (Q]|2]i, we find that their large-scale component is 
described by 



[d, - vV 2 ) u = -(u V)u + (b ■ V)b - (u' • V)u' + (b' ■ V)b' - Vp + f, 
(<9,-?/V 2 )b = Vx(uxb) + Vx(u' xb'), 



(287) 
(288) 



with, in addition, V ■ u = V ■ b = V ■ u' = V ■ b' = 0. The effect of MHD turbulence on u is contained in the 
terms — (u' ■ V)u' + (b' ■ V)b', corresponding to the Reynolds-Maxwell stress-tensor Rjj = -u^u'j - b'.b'.. The effect of 

MHD turbulence on b is contained in the term V x |u' X b'J, where u' x b' is the so-called mean electromotive force 
dSteenbeck et alj,[l966k . 



4.3.2. Reynolds-Maxwell stress-tensor 



The first attempt to estimate the Reynolds-Maxwell stress-tensor using a shell model was made by Geertsema and Achterberd 



(119921) , in the context of a thin, differentially rotating, accretion disk. A vectorial (instead of scalar) shell model was in- 
troduced. Using a local system of cartesian coordinates (x, y, z) in the radial, azimuthal and vertical (normal to the disk 
plane) directions, they introduced three variables per shell for each field , U„ = (£/*, U y n , U„) and B„ = (B*, B y n , B|). 
The corresponding shell equations become 



U„ = W„(a,U,U)-W„(a,B,B)-y^U„+P/. 
B„ = W„(a, U, B) - W„(a, B, U) - ^B„ - ^QB*e y , 



^£lU;e y - 2Q x U„ 



where 



W„(a, X, Y) = k n d [(a„ • X n ^)Y n ^ - A(a n+l ■ X„)Y n+1 ] + lc„C 2 [(a n -i ■ X„)Y„_! - A(a n ■ X n+1 )Y„ +1 ] 



(289) 
(290) 

(291) 



is a vectorial generalization of model ( 1176b , and a„ is a set of arbitrary unit vectors. In Eg. d2891>. P„ is an ana- 
logue of the projection tensor of Eqs. (137H38I I though here it applies to linear terms only ( Geertsema and Achterberj, 
19921) . The Coriolis forces -2Q x U„ are also included. The terms |£2U;fe y and -|fiB ( *e y are the contributions 
from Keplerian differential rotation. Numerical solutions show that the turbulent shear-stress £ [BIB], - U*U)\ can 

n 

supply the strong effective d issipation needed to explain the dynamics of accretion disks as originally suggested by 
Shakura and Sunvaevl dl973l) . 



4.3.3. A subgrid shell model 

A subgrid shell model was introduced bv lFrick et al.1 (120021) for HD convection in a rotating spherical layer heated 
from below, with application to the Earth's core in mind. The large-scale flow u and temperature 9 satisfy 



d,u + (u ■ V)u 

dfi + u- V(0 + Go) 



-Vp + (v + v,)Vu + f, 
(k + k,)V 2 9, 



(292) 
(293) 



where f includes the Coriolis and Archimedean forces, and Go is the temperature profile prescribed throughout the 
layer. In addition, appropriate boundary conditions are applied. The system of Eqs. d292H293b is solved by DNS 
with the resolution given by the grid size. The effect of subgrid turbulence on the scales larger than the grid size will 
be modeled by the turbulent transport coefficients v, and K t . For scales smaller than the grid size the GOY model 
given by Eq. (1185l l for e — 1 - /V 1 (3D turbulence) is solved. In order to provide the correct linkage between the 
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DNS and the shell model, the DNS kinetic energy is calculated in different Fourier shells as introduced in Sec. 12.2.41 
except that here the sequence of shells is geometric. This provides the velocity U„ in the two first shells of the shell 
model, in which no other forcing is applied. The dissipation rate e in the shell model is the total energy dissipated 
per unit of time. It leads to numerical values for v, and k,, estimated as v, * k, ~ 0.1(^e) 1/3 , where l c is the mean 
scale corresponding to the two common shells between the DNS and the shell model. In Fig. [26] the kinetic energy 
is plotted versus the shell number. The two common shells are n — 4 and n = 5 indicated by black squares. Such 
an exam ple demonstrates the feasability of using a shell model as a sub-grid model. Combined with the previous 
example ( Geertsema and Achterbergl 1992 ) it could be used for the creation of a vectorial subgrid MHD shell model. 




Figure 26: Kinetic energy versus shell number, from DNS (filled circles) and shell model (empty squares). The straight line indicates a Kolmogorov 
slope of -2/3. The two common shells correspond to the filled squares. From lFrick et all 120021) . 



4.3.4. The mean electromotive force 

The theory of astrophysical dynamos has been developed mainly using the framework of the mean-field approach 
dKrause and Radlerl [l98C " 



riy sica l dynamos has been developed mainly using the rramework or the mean-held approach 
a lRadlerll2007albh . The latter yields various models for large-scale magne tic fields of celes- 



tial bodies, such as galaxies, stars and planets, reproducing to s ome extent the available observations dRuzmaikin et all 



198 



2003 



Beck et al. ]l996tlB~randenburg and Subramanianl 2005 ). This approach has recently been cha llenged dVishniac et al. 
Cattaneo and Hughesll2009l) and its compatibility with small-scale dynamo action questioned ( Cattaneo and Hughesl 



200 lh. At the heart of th e controversy is (i) how the mean electromotive force u' x b' (the mean e.m.f.), is estimated 
dCourvoisier et al. , 2010l) and (ii) if it is at all relevant to the existence of large-scale magnetic fields. 

Here we assume that the large-scale magnetic field is generated only by the mean e.m.f. (u = 0). The large-scale 
ma gnetic field b is decomposed into its poloidal and toroidal parts, respectively related to the scalar quantities bp and 
b T dMoffattL[l978l) . The latter satisfy 



{d t + nk 2 L )b P = {Vx(u'xb')} p , 
(d.+Tjk^bT = {V x (u' x b')} r , 



(294) 
(295) 



where ki is the wave number of the large-scale magnetic field, and the terms on the right hand side (r.h.s.) are the 
scalar quantities related to the p oloidal and toroidal parts o f the mean e.m.f. The idea of estimating the r.h.s. terms 
using a shell model started with ISokoloff and Frick d2003l). The r.h.s. terms have been deduced e ither using some 
parametrization derive d from the mean-field a pproach dSokoloff and FrickL l2003uFrick et al. . 120061) . or directly from 
the shell model itself dNigro and Veltril 1201 lh . As for the subgrid shell model, the essential point is to provide the 
correct linkage between the large-scale equations and the shell model. 
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In the mean-field approach, the terms on the r.h.s. of Eqs. 
effects, such that 



1294H295t are parametrized with the so-called a and B 



{V x (u' x b')} p * ik L ab T -Bk 2 L b P , {V x (u' x b')} r ~ -ik L ab P -Bk 2 L b T , 
where, using approximations, the a and B parameters are estimated to be 



a — a + a 



a" = -- <u' ■ V x u') , 



I(b'-Vxb'), P*Uu' 2 ), 



(296) 



(297) 



t being some characteristic time of turbulence. The term a", is the usual kinematic tr -effect and acts in favour of 
dynamo action. The term or is related to the feedback effect through the Lorentz forces (iFrisch et all 1 19871) and acts 



against the dynamo. The term B is the so-called turbulent magnetic diffusivity. With a shel l model of L2-typ e (GOY 



or Sabra), and identifying r with the turn-over time, these terms can be written in the form (IFrick et aU 120061) 



(298) 



Sokoloff and Frickl(l2003l) took a h = 0, but included other feedback effects related to Alfven waves. In lNigro and Veltril 
(EoT l|), the terms on the r.h.s. of Eqs. ( I294H295I ) are estimated as ik L (U*B„ - U„B*), where U„ and B„ are obtained 
from a shell model. 

Whatever the expression of the mean e.m.f. in terms of U„ and B„, when considering the full model composed of 
Eqs. d294H295l ) plus the shell model equations, at least the total energy should be conserved in the absence of viscosity 
and diffusivity (v = r] = 0). This leads necessarily to additional terms in the shell model equations. For details we 
refer the interested rea der to the previous ly mentioned papers. 

In their shell model Frick et al. (2006)' used a random forcing on two scales to control both the injection of energy 
and kinetic helicity. The rate of energy injection was constant e = 1 and the rate of kinetic helicity injection was 
varied £ = 0; 0.04; 0.16. In Fig.[27]the kinetic and magnetic energies E u and E B are plotted versus time, together 
with the large-scale magnetic energy E b - \(b 2 p + bj). For v = rj = 10~ 6 the results are given from top to bottom for 
increasing values of £, and from left to right for decreasing values of k^. 

The first observation is that the time span corresponding to the kinematic growth is always roughly the same for 
small-scale magnetic energy (green curve), whereas for the large-scale magnetic energy (black curve) it decreases 
with ki and f . This is directly related to the growthrate of the large-scale magnetic energy which can be estimated as 
Ekin ~ aki -(J3 + rj)k 2 L . In the kinematic range, a « a" cc g. Then for small ki, Eki n ~ gki, implying that Ek; n increases 
with ( and k^. The second observation is that the level of saturation of the large-scale magnetic energy (black curve) 
increases with both £ and k L l . This is related to (i) the definition of a" given by Eq. ( |297t , which increases with ( and 
(ii) to the magneti c dissipation which is proportional to k\. 

Taking a b = O. lFrick et al. I d2006l) calculated the cross-correlation between a" and the large-scale magnetic energy 
E b in the saturated state. They found a systematic delay between a" -quenching due to the growth of E b , and the 
growth of E b due to a". Such dynamics excludes any simple algebraic relation between a" and E b , in contrast to what 
is usually assumed in the mean-field approach. 

In Fig. [28] the mean values of the three energies E u , E B and E b versus Pm are plotted for a given value of the 
viscosity v = 10~ 6 . For Pm < 10~ 4 , the small-scale dynamo works poorly and even stops for Pm < 10~ 7 because the 
magnetic Reynolds number becomes too low. Thus for Pm <K 10~ 4 , \a b \ <sc 1, implying that the large-scale magnetic 
field is generated by the a"-effect alone. Note that the level of large-scale magnetic energy is even greater than that of 
the kinetic energy. Now as v is fixed, increasing Pm corresponds to raising Rm thus favoring the small-scale dynamo. 
If so then |ar fo j also increases, implying quenching of the large-s cale dynamo as depicted by the negative slope of 



E b versus Pm in Fig. [28] Similar results were obtained by DNS dPontv and Plunianl 1201 lb . The former results for 
Pm < 10 4 would presumably drastically change if Pm rather than v was kept fixed. Indeed the reason why the 
small-scale dynamo shrinks for Pm < 10~ 4 is that Rm becomes low. Increasing Rm, while keeping Pm fixed at a low 
value, might lead to a stronger small-scale dynamo and to quenching of the large-scale dynamo due to an increase of 



For a given viscosity v = 10 . iNigro and Veltril (1201 11) studied the dynamics of their multi-scale dynamo intro- 
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Figure 27: Time series for E b (thick black curve), E B (green curve) and E u (red curve). From top to botto m, the inject i on rat e of kinetic helicity is 
increased, f = 0; 0.04; 0.16. From left to right, the scale separation is increased, kt = 1 /2; 1/8; 1 /32. From lFrick et alj 120061) . 



duced above. In particular, they found a hysteris cycle represented in Fig.[29]for both ratios E b /E u and E B /E u versus 
Rm. On increasing Rm, the dynamo starts for Rm as 65. On decreasing Rm, the dyn amo remains for Rm < 60 . Such 



sub critical dynamos are a lso observed in DNS, e.g. for rotating convective dynamos (iMorin and Dorm vL 12009) 



Step anov et al. d2006h . with galactic disks in mind, elaborated a multi-scale dynamo model in which the toroidal 



part of the large-scale magnetic field is generated by differential rotation rather than by means of mean e.m.f.. In 
the mean-field terminology the dynamo is denoted as an aw-dynamo, in contrast to the previous models denoted as 
^-dynamos. The large-scale magnetic field equations are 

d,b P = ak L b T +pk L I d 1 z bp, (299) 
d,b T = -Gd-bp+pk 2 L dlb T (300) 

where bp and bp not only depend on time but also on the coordinate z perpendicular to the plane of the galactic 
disk, with -1 < z < 1. Appropriate boundary conditions in z are applied to bp and bp and their first z-derivatives. 
The parameter G is related to the intensity of the differential rotation. The parameter a also depends o n z as a 



(a" + or) sinfe). The linkage between the shell model and the large-scale model is similar to that in iFrick et al 



(12006b . In contrast to the previous models this model has to be integrated in z. In Fig|30]the z-profile of the poloidal 
and toroidal magnetic field components are plotted versus time. The poloidal component exhibits a much stronger 
dependence on z (smaller scales) and varies faster than the toroidal component. The toroidal component exhibits one 
reversal. 



4.3.5. Reversals of large-scale magnetic field 

The system of Eqs.([T]|2]i is invariant on changing b to -b, allowing for opposite magnetic field polarities with 
the same velocity field u. Magnetic reversals are indeed observed, e.g. the Earth geomagnetic dipole fluctuates 
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Figure 28: E u (white squares), E B (black diamonds) and E (crosses) versus Pm,forv= I0~ 6 , f = 0.08 and k L = 1/16. From, Fric ketafl 120061) . 



and reverse s chaotically ( Merrill and McEl hinnvl 119831) . while the s olar magnetic field re verses with a rather stable 
periodicity ( Weiss and Thompson! 2009 ). The experimental results of Ravelet et al. (2008) showing chaotic magnetic 
reversals have recen tly emphasized the impor t ance of small - scale turbulence in triggerring the reversals of the large- 
scale magnetic field (IPetrelis and FauveL 1200 8: Petreli s et al 1 120091) . It is then rather natural to expect chaotic reversals 
in multi-scale shell models. This is visible in Fig. [27] where each time the black curve, corresponding to E b , crosses 
the horizontal axis, b is changed to —b. 



Ryan and S arson! d2007l 1201 ll) considered a large-scale <m>-dynamo model in which a is calculated from the 



3D HD GOY model (no small-scale magnetic field, no back-reaction from the large-scale magnetic field onto the 
turbulence). They find that the reversals fit a log-normal distribution as do the paleomagnetic data, stressing the 
importance of multiplicative noise in the underlying geodynamo. 

Co nsidering the MHD Sabra model given by Eqs. (1237H2381 I. iBenzi and Pintonl d2010h and iNigro and Carbonel 
d2010l) added a term to the magnetic field equation for one of the largest shells, say n = 1, such that the new equation 
for B\ becomes 

d,Bi = Wi(U,B)-Wi(B,U)~ M(Bi). (301) 

They chose M(Bi) oc B\ or M(B\) oc B\ ( Benzi and Pintonl 2010l) . or a combination of both ( Nigro and Carbone , 
2010h . Of course energy conservation is thus violated, but this was the simplest way of breaking the symmetry of the 
system and have reversals for B\ with p eriods of constant sign. T hey found that the mean period between reve rsals 
increases with the magnetic diffusivity ( Benzi and Pinton , 2010l) and the viscosity ( Nigro and Carbone , 2010t) . In 
Fig-Elan example is given for different values of viscosity. From top to bottom the viscosity is decreased by a factor 
10, while the mean period between two reversals clearly becomes shorter. 

4.4. Alfven waves and rotation 

Assuming that MHD turbulence is embedded in a homogeneous magnetic field bo, and in a rotating frame with 
angular velocity £1, Eqs. ([T]|2]i are changed into 



(d, - vV 2 ) u = -(u ■ V)u + (b ■ V)b + (b • V)b - 2Q x u - Vp + f, 
[d, - 77V 2 ) b = -(u V)b + (b ■ V)u + (b ■ V)u, 



V • u = 0, 

V • b = 0, 



(302) 
(303) 



As explained in Sec. 12.1.11 Alfven waves may propagate along bo in both directions provided that dissipation is not 
too str ong. This later statement is actually what makes experimental evidence difficult to obtain, especially in a liquid 
metal (lAlboussiere et al.L 1201 11) . For a sufficiently strong bo or il the MHD turbulence becomes anisotropic, making 
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Figure 29: Ratio (top) and E B /E U (bottom) versus Rm for v = 10~ 5 . Adapted from lNigro and Veltril J201lh . 



the use of isotropic 3D shell models questionable. In order to account for such anisotropy, other shell models have 
been elaborated as presented below. Mainly three lines of research have been followed so far. 



Two shell models of Alfenic turbulence have been elaborated by ICarbone and Veltril (11989 , 1990) within an 
isotropic and anisotropic framework. In the latter case they introduced an angular dependence with respect to the 
direction of anisotropy (the direction of bo). The model coefficients were estimated using a standard statistical 
approach (Direct interaction approximation, Markovian approximation and random phase hypothesis). The 
interested reader should refer to the above mentioned papers. 



Another line of developmen t (Sec . 14.4.1b co nsists in keeping an isotropic MHD shell mode l and simply adding 
a term correspond i ng to b p ( Biskamd," 1994 Hattori et all 200 lb . or SI for HD turbulence ( Hattori et al. , 2004 ; 
Chakrabortv et al. , 2010t) . or both jPlunian and Stepanovi l2010b . Such a model is, of course, well suited to 



isotropic phenomenology like IK for Alfen waves. 

The usual way of dealing in the presence of an applied magnetic field bo is to write the MHD equations in the 
plane perpendicular to the direction of bo, leading to the so-called reduced MHD equations (RMH D). The first 
RMH D shell model was developed in the context of intermittent heat ing in solar coronal loops (|Nigro et al 
20041). 'ts statistical properties compare successfully with observation (Nig ro et all l2005t iBuchlin et all 120051 



Buchlin and Vellil l2007t iBuchlinl 12007b . In the context of the Alfvenic solar wind the RMH D shell model gi ves 
a good description of the transition between weak and strong MHD turbulence (IVerdini and GrappirA 120121) . It 



also provides a mec hanism for the presence of the low freque ncy magnetic spectrum observed inside the sub- 
Alfvenic solar wind ( Verdini et al. i l2009t Verdini et al., 2012a). The model and a few results are summarized in 
Sec. l4A2l 

4.4.1. Isotropic shell models 

A shell model with an externally applied magnetic field and global rotation of intensities bo and £2, can be written 
in the form 



d,U„ = W n (U, U) - W„(B, B) + ik„b B„ + iQt/„ - vkiU„ + F n , 
d,B n = W n (V, B) - W n (B, U) + ik„b U„ - i]klB„. 



(304) 
(305) 



Obviously, taking bo + and Q. + implies the failure of magnetic and cross helicity conservation respectively, in 
agreement with the original equations. In addition to the eddy turn-over time t^L = U u h the introduction of bo and O 
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Figure 30: Isolines of poloidal (to p) and toroidal (bottom) component of the large-scale magnetic field versus z (vertical axis) and time (horizontal 
axis). From lStepanov et al.1 120081) . 

leads to two other time-scales, t\ = I /bo and = Q. 1 , over which energy transfers may occur. Depending on which 
is smallest, f^z,, tA or tn, one expects different types of turbulence characterized by different slopes for the energy 
density spectrum within the inertial range, respectively E(k) oc £~ 5 / 3 (type K), E(k) oc k^ 3 ^ 2 (type A) or E(k) oc k~ 2 
(type R). As fjvx and depend on scale /, an inertial range with several slopes generally occurs. 

In the top-left panel of Fig. [32] we summarize the range of values of parameters bo and O for which the four follow- 
ing regimes are possible, provided Pm = 1, Rotation (R), Rotation-Kolmogorov (RK), Rotation-Kolmogorov-Alfven 
(RKA) and Rotation- Alfven (RA). Each regime is characterized by different inertial ranges with, going from small to 
large k, slope k~ 2 for (R), slopes k~ 2 and k~ 5 ^ for (RK), slopes k~ 2 , k~ 5 ^ and k~ 3 l 2 for (RKA), and slopes k~ 2 and Ar 3 ^ 2 
for (RA). In addition, two other regimes are also possible for Q. — 0, Kolmogorov (K) and Kolmogorov-Alfven (KA) 
depending on whether bo is weaker than (rje) 1 ^ 4 or not. In the other panels of Fig. [32]sets of results are shown for the 
cases (RA), (RK) and (KA) with v = 10~ 7 and Pm = 1, using the helical LI shell model given by Eq. (12641 i. We find 
excellent agreement between the results obtained with the shell model and isotropic phenomenological predictions. 

For Pm «: 1, the characterization of the different regimes is more complex, but still tractable analytically. How- 
ever, the results obtained with the shell model clearly show that the distinction between the different regimes is in most 
cases no longer possible. This is because the magnetic dissipation scale for Pm < 1 is larger than the Kolmogorov 
scale implying that the spectral slopes are not c lear e nough. 

Finally, note that in Plunian and Stepanov ( 2010h it was necessary to introduce correlation times for bo and £1 in 
order to avoid spurious supercorrelation between U„ and B„, as previously explained (Sec. 14.1. Ti l. This is presumably 
related to the number of degrees of freedom in the helical LI -model ( 1264b . For HD turbulence with a constant Q we 
verified that such a problem disappears using the helical L2-model given by Eq. ( 1269l >. 

4.4.2. RMHD shell models 

The Reduced MHD equations accounting for anisotropy due to a strong ext ernally applied magnetic field bo, are 
obtained by restricting the initial MHD Eq. (|4]i to field-perpendicular equations (iStraussl 1 19761) 



3,z± + bod x z* + (zl ■ VJz* + V ± p = r + Viz ± + rV 2 ^, V ± ■ = 



(306) 



where x is the coordinate along the direction of bo, is the projection of the Elssasser variables onto the plane 
pe rpendicular to bp, an d r* = + rj). 

iNigro et al.l (120041) introduced a RMHD shell model for Eq. (I306l l keeping the x-dependency along the direction 
of bo (with appropriate boundary conditions) while using a 2D shell model in the perpendicular direction to account 
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Figure 31: Large-scale magnetic field versus time. From top to bottom the viscosity is decreased. Adapted from[Niero and Carbone J2010t> - 



for 2D MHD turbulence. The spectral space is thus divided into cylindrical shells, each shell n being characterized by 
k„ < |kj < k n +i, where k ± is the component of the wave number perpendicular to bo- The shell variables U„ and B n 
depend on both time and x. Such an hybrid shell model written in Elssasser variables Z* = U„ + B n , is given by 



(d t + bod x )2% = W n {Zl,Z±) - kl (r + Z± + r~Z*) ., 



(307) 



Nigro et al.l d2004 used the MHD GOY model ( fT92l with A = 2 and e = 5/4 to obtain 



~ /13 11 19 11 19 13 \* 

w„(z:,z^ = +ik n y-zt, 2 z: +l + -zj^, - —z^^ - ^ +1 zt, + ^K-^u - ^ z «-^-ij ■ vw) 

To model a coronal loop, such as the one represented in the left-panel of Fig. [33] kinetic energy is injected at 
one foot x — of the coronal loop, with no motion in the other foot x — I. At x = the motion is non-zero only 
in the first three shells (« = 0,1,2) and has a Gaussian time distribution. With appropriate values for the different 
dimensional pa rameters (motion a t x = 0, length and radius of the coronal loop, Alfven wave velocity), and for 



io- 



Nigro et al. (2004) estimated the time evolution of the total energy, net incoming energy flux and 
dissipated power (right panel in Fig. [33l l. This study strongly supports the idea that coronal nanofiares are due to 
intermittent events produced by Alfven wave turbulence. About 60% of the energy entering the system is in fact 
dissipated. The dissipation compares well with the energy involved in nanofiares and, from Fig. [33] clearly displays 
intermittency with statistics comparable to those of nanofiare emissions. Introducing a more complex model with 
several layers in x. lVerdini et al. (2012b) studied the combined effects of turbulence and energy leakage on the coronal 
heating. 



The same RMHD shell model given by Eqs. ( 130711308b was also used bv lVerdini and Grappinl (120121) to investigate 



the transition from weak to strong turbulence, depending whether t& <s fjvx, or « t^i, where oc (k\\bo) is 
the characteristic time responsible for Alfven waves propagation, and t^i K (k±b(k ± ))~ l is the eddy turn-over time 
(see Sec. 12.2.2b . In order to prescribe the ratio f^/fiVL , they added a forcing term /* (x, t) in Eq. ( 1307b leading to an 
estimation of and fjvx at the forcing scale. Typical perpendicular spectra are shown in the left-panel of Fig. l34lfor 
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Figure 32: Top-left panel: all possible non-linear regimes in the map (bg, fl). Other panels: normalized energy density spectrum versus normalized 
wave number. They correspond to reg imes RA (top-right), RK (bo ttom-left) and KA (bottom-right). The different curves correspond to different 
values of bo, il or both. Adapted from lPlunian and StepanovlfeOlCl) . 



tA/tNL ~ 1 an d ~ 1/32. They found that (a) only one slope k' 5 ^ 3 is obtained for strong turbulence, whereas (b) two 
slopes Ic 1 and lc 5 ^ are found for weak turbulence. In the notation introduced in Sec. 14.4.11 the weak regime is of 
(AK) type, in striking opposition with the (KA) is otropic case. 

In the context of solar wind MHD turbulence, IVerdini et al. (2009) introduced a modified RMHD shell model in 
the form 



dZ± dZ± 1 (ldp 

—!L + (U±b Q )-^--(U + bo)[--f- 
at ox 4 \p dx 



1 (I dp IdR 

+ -(U + bo)[ — - + 2 

4 \p dx R dx 



Zl = W n (Z*,ZZ)-% (r + Z± + r-7%) , (309) 



where U, p and R stand for the background solar wind, the fluid density and the radius of the magnetic flux tube 
respectively. The profile of the latter quantities as a function of x was assigned along with bo. The modulus and 
correlation time of Z+ was prescribed at x — and only in the three first shells. This corresponds to the in jection 
of high frequency Alfven waves at the base of the chromosphere. By solving Eq. d309t IVerdini et al.l d2012al) found 
a double power law for the magnetic frequency spectrum that is reproduced in the right panel of Fig. [34] At low 
frequency it comp ares well with t he 1 // magnetic spectrum of the sub-Alfvenic solar wind represented in Fig. [3] For 
higher frequencies Verdini et al. ( 2012al) argue that the -2 slope that they find in their model is masked in the solar 
wind by the more energetic perpendicular spectrum which has a -5/3 slope. 



4.5. Hall-effect 

In strongly magnetized conducting media, e.g. in weakly ionized accretion disks, white dwarfs and neutron stars, 
the Hall drift of the magnetic field can operate on a much shorter time-scale than the other transport processes. Though 
non-dissipative, the Hall current redis tributes the magnetic energy from large to small scales and thence enhances the 
dissipation rate of the magnetic field ( Urpin and Shalvbkov , 19991) . The rate of magnetic field decay is important to 



62 




\ J "I 




<b) 


- 1 Ulylill. 


(0 



S 10 IS 20 25 30 35 40 45 SO 55 80 65 70 75 SO 
Time (hours) 



Figure 33: Lef t panel: Coronal loop in the direction of bo- The 2D shell models in planes prependicular to bo are p iled up along bp. A dapted from 
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Figure 34: Left panel: Reduced perpendicular energy spectra normaliz ed by k (solid lines). Th e do t-dashed line is the k , 2 scaling. Right panel: 
Reduced frequency spectra for E"' b at r wp (solid lines). Adapted from[Verdini and Grappin 1 2012 ) and Verdini et al. 1 20123). 



understanding the history of such astrophysical objects. The Hall effect is probably most pronounced in young neutron 
stars, called magnetars because of their strong magnetic field of up to 10 15 G. 

In the absence of motion and am bipolar diffusion, the evolution of t he magnetic field due to the Hall effect can be 
written in its dimensionless form as (IGoldreich and Reiseneggen . 119921) 



^ = R H 'V 2 b - V x [(V x b) x b] 



(310) 



where Rh = (ebot e )/(m* e c), bo the characteristic magnetic field strength at the largest scale, e the elementary charge, c 
the speed of light, m* the effective mass of an electron and t e the electron relaxation time. 

The last term of Eq. (131 Oi l describes the advection of the magnetic field by the Hall drift and is expected to produce 
a non- linear cascade from large to small magnetic scales. Using the assumption of strong turbulence, Vainshteinl 
(119731) derived a kr 1 ^ scaling law for the magnetic energy density spectrum. However, considering weakly interacting 



waves 



Goldreich and Reiseneggeri i 19921) found an alternative k 2 scaling law. In both cases the spectra are steeper 



than for Kolmogorov turbulence. 
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The possibility of a cascade is of course appealing for shell models. In addition, for Re — > °° Eq. (13 1 Ob gives rise 
to two quadratic invariants, the energy E b and the magnetic helicity H b , which are sufficient for the derivation of a 
shell model including the Hall effect. The difference with MHD shell models is that the non-linear term is now of the 
for m k„Q„(X) 



Frick et al.1 d2003l) introduced a Hall shell model in the form 

-12 i( e 1 — e \ 
(d t + R H k n )B„ = k~yB n+ 2B„ + i - -^B n+ \B„^\ B„_iB„_2 J + F„ , 



with the two quadratic invariants 



(311) 



(312) 



In Eq. d31 U F„ is an externally applied electromotive force and B n is real. Indeed the real and imaginary parts of 
complex variables B„ would not couple, as can be seen from Eq. ( 13101 ) written in Fourier space. We take s = 1 - A in 
order to impose the conservation of magnetic helicity 



(313) 



Solutions of Eq. ( 131 11 1 for F„ = (free-decaying Hall-turbulence) show that the initial magnetic energy concentrated 
in the large scales tri ggers the cascade process. This results in a magnetic inertial range of the form E B {k n ) ~ k„ A ^, 
corresponding to the IVainshteinl (I1973h ener gy density E b (k) ~ k~ 7 ^. However, the energy spectrum steepens rapidly 
in time, unless stochastic forcing is applied (IFrick et all 120071) . 

On representing the magn etic field as the sum of poloidal bp and toroidal bj components, Eq. d3 10b takes the form 
dUrpin and Shalvbkovi \l99% 



ReVbp- V x [(V x b T ) x bp] 



<9b P 

Jj^ = R H 1 V 2 bT-Vx[(Vxb T )xbT + (Vxbp)xbp] 



(314) 
(315) 



From Eqs. ( 1314113151 . we clearly see that the two components bp and bj are coupled. However, coupling is not 
symmetric. An initial poloidal configuration can generate a toroidal component with the term -V x [(V x bp) x bp]. 
On the other hand, a poloidal component cannot be generated if the initial magnetic configuration is purely toroidal 
unless some transient dynamo type instability occurs for bp. Note that, in t he presence of a background m agnetic 
field, such a mechanism is at the heart of the so-called Hall drift instability (IRheinhardt and Geppern , 120021) . In the 
context of neutron stars, this instability leads to non-local energy transfer to small scales and then to enhanced crustal 
fiel d dissipatio n . 



Frick et al.l d2003l) derived the following shell model for both poloidal P and toroidal T components of the magnetic 



field 



d t P„ + R^k^P „ — kf t \ P n +2T n +i - 



d t T n +R^k 2 n T n = k 2 n \P n+2 Pn + \ 



A+l 1 

, 2 Pn+U n -\ + — ^P n -{T n -2 \ + k n I T„ + 2Pn+l 

A+l 1 

' , 2 + -pPn-\Pn-2 I + k n \ T n+ 2T,,+ \ 



with magnetic energy and magnetic helicity defined by 

E^^iPl + T?,), H B = J]k^P n T„. 



— "j-j— TWl^n-l + ^r„_iJV2j,(316) 

- ^-r„ +1 r n _! + -^r„_,7v 2 W317) 



(318) 



As in helical models a remarkable advantage of definition (131 81 > is that, contrary to (1312l i. magnetic helicity and 
magnetic energy are not shell-by-shell correlated. Another advantage is that it is possible to describe the energy 
exchange between both poloidal and toroidal components of the magnetic field. In particular, it was shown that a 
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transient poloidal field component can develop from an i nitial purely toroid al configuration ( IFrick et all 120031) . With 
again stochastic forcing incorporated into Eqs. (l3 16B3 17t jFrick et al.l ( 12007b found a stable inertial range characterized 

by the IVainshteinld 1973b energy spectrum slope E b (k) ~ Ic 11 '-' and H b ( k) ~ kr l E b (k). 

With the application to the solar wind in mind, Hori et al. ( 2005 ) and Galtier and Buchlin] ( 2007 ) included the Hall 
effect in the incompressible MHD Eqs. (Q]|2]i 



(<9, -vV 2 )u = -(u ■ V)u + (b • V)b - Vp, V ■ u = 0, 

(<9,-?7V 2 )b = -(u- V)b + (b- V)u-R H Vx [(Vxb)xb], V b = 0. 



(319) 
(320) 



In the limit (v, 77) — > (0,0), the system of Eqs. d3 19113201 ) contains three quadratic invariants, the total energy E and 
magnetic helicity H" defined in Eqs.((5]|6jl, and a third (new) quantity 



"'-I 

Jv 



(a + R H u) ■ (b + R H V x u)dV, 



(321) 



called the ion helicity. Note that H' = H b + 2R H H C + R H 2 H". In the limit R H -> 0, H' reduces to magnetic helicity, 
while (H' - H b )/2Rii, also a conserved quantity, reduces to cross helicity H°. 



In the spirit of the complex G OY shell model, the Hall MHD shell model for A = 2, takes the form (IHori et al 



2005; Galtier and Buchlin, 2007) 



(d, + vkl)U„ = ik n [(U n+1 U n+2 - B n+l B n+2 ) - i(t/„-it/„ + i - B„-iB„ + i) - \(U n - 2 U n -x - B,^ 2 B„^)f+F„. 



(d, + T]kl)B n 



ikn 
6 



U n +\B n+ 2 - B„ +1 U„+2 + U n -\B n+ \ - B„_it/„ + i + t/„_2-B n -i _ B„-2U„- 



n • 2 [ ^ 1 * 

+ (-\)"iRYik n \B n+ 2B„+\ - -B n +\B n -i - -B„-iB„-2} ■ 
In the limit (v, if) — > (0, 0) the following quantities are conserved 

£B= 2~X |Bn|2 ' = \Yi { ~ l) " k » l \ B "\ 2 > H 1 = H B + 2R H H C + R H 2 H L 

n n 

with 



(322) 



(323) 



(324) 



(325) 



Note that in Eq. ( 1323b the term accounting for the Hall drift differs from Eq. ( 131 11 1 because here the variables are 
co mplex. 



Galtier and Buc hlin (2007J) and lHori and Miu ra (2008) found that the large-scale magnetic energy density follows 



a k~ 5 ^ spectrum that steepens to k~ 1/3 at scales smaller than Rh if the magnetic energy overtakes the kinetic energy 
or to k~ n ^ in the inverse case. This might explain why the magnetic spectrum of the solar wind steepens at high 
frequencies (/ > \Hz). 



5. Summary and outlook 

Table [3] summarizes the MHD shell models that have been discussed in the present review. The L2 GOY model 
given by Eqs. ( 1192H194t is without doubt the most often used, in its 2D and 3D forms. It gives excellent results in 
terms of spectra, intermittency, and energy transfer and flux. However, we saw that such a model can be misleading 
when dealing with kinetic helicity in 3D HD turbulence. This is also true for its cousins the L2 and N2 models given 
respectively by Eqs. ( 123711238b and Eqs. ( 123411236b . On the other hand helical models do not suffer from this drawback. 
So far only the HI local model given by Eq. ( 1264b has been used for MHD turbulence. The other helical models given 
by Eq. ( 1268b or Eq. ( 1269b are given here for completeness, but have not been tested yet. 
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As shown in Appendix, the models depicted in Table [3] represent only a subset of all possible models that satisfy 
total energy and cross helicity conservation. In addition any combination of models will again satisfy both conserva- 
tion laws. In Table[3]the models which depend on one free parameter are the result of the combination of two other 
models. Regarding the non-local parameter in Table[3] it controls the degree of non-locality in non-local models. We 
saw that this parameter can be estimated with the help of a hierarchical approach or from phenomenological arguments 
(Sec. 03). 

In addition to automatically satisfying the conservation of total energy and cross helicity, the general formulation 
given in Eqs. ( 114511147b has the advantage of offering a mathematical framework for the definition of energy transfer 
and flux. This general formulation has also been used to define transfer and flux of kinetic helicity in HD turbulence. 
No doubt it can also be applied to transfer and flux of cross and magnetic helicity in 3D MHD turbulence. 

The dimension of MHD turbulence, 2D or 3D, can be changed imposing either the square of magnetic potential 
or the magnetic helicity as the third quadratic invariant. Though this is possible to achieve for any models, except LI 
for which magnetic helicity cannot be defined, in Table|3]we denote by 2D, 3D or both the corresponding version that 
has been developed so far. 

Finally modern shell models always use complex instead of real variables mainly because it increases the degree 
of freedom of the system, leading to more realistic dynamics. 

The applications presented in Sec.|4]correspond to a large panel of situations requiring high values for the kinetic 
and magnetic Reynolds numbers: small-scale and multi-scale dynamos, free-decaying MHD turbulence, Alfven wave 
turbulence and Hall-effect turbulence. In order to deal with anisotropy hybrid models have been developed by mixing 
2D MHD turbulence in a plane perpendicular to the direction of anisotropy with direct integration in the parallel 
direction. A few attempts have been made to develop subgrid shell models, which are indeed an interesting and 
promising application of shell models. 

Shell model simulations require the same time-stepping as direct numerical simulations. However, there is a 
significant gain in using shell models, essentially because the number of grid points is much lower than in DNS, 
and partial derivatives are absent. Shell models are user-friendly tools that guide our intuition in realistic parameter 
regimes still inaccessible to DNS. 



Acknowledgements 

This collaboration benefited from the International Research Group Program supported by Perm region Govern- 
ment. R.S. acknowledges the grant YD-447 1.20 11.1 from the Council of the President of the Russian Federation 
and the RFBR grant 1 1-0 1-9603 1-ural. PF acknowledges the RFBR grant 1 1-01-00423. We warmly acknowledge J. 
Torbet for great improvement of the writing. 



Appendix A. Ll-models 

The complex Ll-models have the form 

+i 

w„(x, y) = k J] 4 x "+< y "+j + b 7jK+?n+j + cfjX n+i Y; 1+j + d%r n+i r n+J , , (A.i) 

i,j=-l,\i-j\<l 

involving 28 complex coefficients. The number of complex coefficients reduces to 1 1 on applying the property (1147l >. 
Then the general shape of complex Ll-models becomes 

w„(x, Y) = k„ [ A l (x* n _ 1 r n _ l - AT n r n+ d + MKK-i - 
+ a 3 (x„_ 1 f,;_ 1 - Ax n r n+l ) + a 4 (x„Ci - 

+ A=,X* n X Y n -i - AA* 5 X„Y„+i + A(,X*Y„-\ - AAlX n+ iY n+ \ 
+ AjX„-[Y„-\ - AAjX*Y n+ i + A%X n Y n -\ - AA\X* n+x Y n +\ 

+ (A 9 X„_i - A* g X*_ l +A W X„ -A* m r n +A n X n+l - A* n T n+l )Y n ] . (A.2) 
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Model 


Equations 


Variables 


Conservative quantities 


Comments 


Reference 


LI 


(I176H178D 


real 

U„, B„ 


E = { %(U 2 n + Bi) 


One free 
parameter 


Gloasuen etal. (1985) 






LI 


(I180H182D 


complex 

U n ,B n 


E = 1 Z(l^l 2 + 

// c = i + u;B n ) 


One free 
parameter 


Biskamo ( 19941 






L2 
GOY 


<U92H194D 


complex 

U„, B„ 


E = \ £(|{/„| 2 + |B„| 2 ) 
# c = \ UU„B: + U*B n ) 

h b = lz(-iyk- l \B n \ 2 

A=fZk- 2 \B„\ 2 


=> 3D MHD 
=> 2D MHD 


Brandenburg et al. (1996) 

Basuetal. (1998) 

Frick and Sokoloff( 1998) 


L2 
Sabra 


(I237H238I) 


complex 
U„, B„ 


E = | Sdf/J 2 + \B„\ 2 ) 
H c = i ZiUnBl + U* n B n ) 

H B = ±Z(-l)"k7, l \Bn\ 2 


3D MHD 


Plunian and Stepanov (2007) 


N2 
Hierar- 
chical 


(12141). 
d220H221b 


real 

U„, B„ 


E = £ Z(u 2 + Bi) 

H c = 2 f/„B„ 


2D MHD 


Frik (1984) 


N2 


(I234H236I) 


complex 
t/„, B„ 


e = \ 2(if/„i 2 + m 2 ) 

H c = i £(U„B*„ + U* n B„) 
H B = lZ(-m; l \B n \ 2 


One non-local 
parameter 
3D MHD 


Plunian and Stepanov (2007) 


HI 


(12641) 


complex 

U n , B„ 


E=\Z(\Un\ z + \B„\ z ) 

h c = \ Z(u n B; + u;b„) 

H B = ±Zk- n HB* n 2 -Bl) 


One free 
parameter 
3D MHD 


Mizeva et al. (2009) 


H2 


(12681) 
(12691) 


complex 


E = \ + \U n \ 2 + |B,tl 2 + IB"] 2 ) 
H c = i 2(C/„ + Br + C/,7B,r + c.c.) 


3D MHD 


This paper 



Table 3: Summary of main MHD shell models. 



where the A,- are complex parameters. After our definition of LI -models given in Sec. l3.1.2l shell « cannot interact with 
its elf only, implying A io = 0. A similar general shape for complex LI -models has been introduced in the seminal paper 
bv lGloaguen et al.l dl985l) but using Elsasser variables Z ± . The Liouville's t heorem yields to additional constraints on 
the possible choice for the A,. The model given by Eq. d 1 80b investigated bv lBiskampI ( 1994 ) corresponds with taking 
A; = 0fori = 3,--- ,11. 



Appendix B. L2-models 

Any complex L2-model has the form 

+2 



W n (X, Y) = k n Yj a?jX„ + iY„ + j + bfjX: +i Y n+j + cfjX^Y^j + <pC + ,T„V (B.l) 

U=-2, 
|i+j|=3 or /=-/'=+ 1 

A total of 24 complex coefficients are involved. On applying the property ( 11471 ) the number of complex coefficients is 
reduced to 12 and the general shape of complex L2-models becomes 

w„(x, Y) = k„ [ c,{x* n _ 2 r n _ x - ax:_ x y: +i) + CiiK-iK-i - ^ 2 K + iK +2 ) + W +1 Ci - ^ +2 Ci) 

+ C\(X n -2^* n -\ - ^n-lT* +1 ) + CiX* n _ x Y n -2 - A 2 ClX n+ \Y„ + 2 + C(,X* +l Y n -i - AC^Xn^Yn+l 

+ C-iX* n _ 2 Y n -\ - ACjX„-\Y n+ i + C$(X„-\Y*_ 2 - A 2 X„ + iY* +2 ) + C<)X„ + iY„-\ - AC\X* n+7 Y n +\ 

+ CmX„-2Y„-i - AC^X^Yn+i + C\\X n -\Y n -2 - A 2 C\^X* n+ ^Y n+ 2 + Cn{X n+ \Y* n _ x - AX n+ 2Y* +l ) ] , 

(B.2) 
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where the C, are again complex parameters. Liouville's theorem is automatically satisfied as W„(X, Y) does not 
depend on X n , X*, Y„ nor Y*. The GOY model is obtained by setting all C, coefficients to zero, except Ci, C2 and C3. 
The Sabra model is obtained by setting all C, coefficients to zero, except C\\,C\2 and C13. 



Appendix C. Nl-models 

The Nl-models are obtained from the Ll-model ( 1A.2I I. including the non-local interactions. Their general shape 



is 



w n (x, y) = k n ^ t Altec, - ^c_ m+1 Ci) + A ?(^C-iC m - ^; +M -i c m ) + - ^ +m c m ) 

m>l 

+ Ag(-2f n _ m Y^_2 — /L2f„_ m+ i Y n+l ) + A^(X n -\Y n _ m — A m X n+m -iY n+rr ^) + A4{X n Y n _ m — A m X n+m Y n+m ) 
+ ^S^n-m^n-i ~ AA^ X n - m+ \Y n +\ + A$X n _ j Y n — m — A Ag X n+m -iY n+m + A( i X n Y n - m — A A (> X„ +m Y n+m 
+ A^X„- m Y„-i — AA^ X n m+l Y n+ i + AjX n -iY n - m — A'"A^j X n+m _^Y n+m + A$X„Y„- m — A ln AgX n+m Y n+m 
+ (A 9 X n - m+2 -A;X*_ m+2 )Y n ]. (C.l) 



where the A, and Ai are complex parameters depending on m. The non-local version of the model ( 1176b bv lGloaguen et al 



(1985) becomes 



W„(X,Y) — k„ ^C\(X„-,„Y„-\ AX n - m+ iY n +\) + C\(X n -\Y n - m — A m X n+m -\Y n+m ) + C2(X„Y n - m — A'"X„ +m Y„ +m )^ , 

(C2) 

where Cj, C\ and C2 are real parameters, and X and Y are real variables. 



Appendix D. N2-models 

The N2-models are obtained from the L2-model ( IB. 21 ). including the non-local interactions. Their general shape 

is 

w n (x, Y) = k n J] [ Ci(x,;-,„-iCi - *K- m Y* n+l ) + ^(i,;,^, - * m+1 K +m K+ m+ + c 3 (x ; ; +1 c„ - i"'x; +m+1 y,; +m ) 

m>l 

, z-" 1 /y V* ]y V*\_i_/'~ 1 V* V im+l^v v _i_ f~* V* V 3 m /"**V V 

t <-4V A n-m-M „_i — / i A n-iii'j + lJ "•" ^5 a ,i-1 j m-1 — A <-5 A n+m-< n+m+l ^-6 A n+ \ I n-m ~ A ^^n+m+X 1 n+m 

+ CiX n _ m _ J F„_ 1 — AC-jX n - m Y n+ \ + Cg (X„_ 1 y„_ m _ 1 — A m+l X n+m Y n+m+l ) + C9X„ + iT„_ m — A m CqX n+m+l Y„ +m 
+ CioX n - m -iY„-\ — AC^X n _ m Y„ + i + C \\X n -\Y n - m -\ — A m+i C\yX n+m Y n+m +\ + C\2(X n +\Y n _ m — A m X n+m+ iY n+m ) ] , 

(D.l) 

where the C, are again complex parameters depending on m. Liouville's theorem is again automatically satisfied as 
W„(X, Y) does not depend on X„, X*, Y„ nor Y*. Non-local version of GOY and Sabra models are obtained by setting 
all C, coefficients to zero, except C\ , C2 and C3 for GOY, and Cn, Cyi and Co for Sabra. 



Appendix E. Numerical aspects 

Computational gain using a shell model 

The computational gain using a shell model rather than a DNS can be estimated. The cost of a simulation is 
proportional to MN, where M is the number of time steps and N the number of grid points. Firstly we consider HD 
turbulence. 

In each direction the number of grid points can be estimated as the ratio between the scale at which energy is 
injected and the scale at which energy is dissipated. Defining the Reynolds number at the forcing scale by Re = 
IfUi f /v, with, from Eqs. ([20b and (l22l , u\ F oc (e/p) 1 / 3 and l v oc e~ 1,,4 y 3 / 4 , we find 

W/ v *Re 3/4 . (E.l) 
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For a DNS in 3D this leads to 

<! S *(W/v) 3 -Re 9/4 . (E.2) 

For an isotropic shell model not only is there just one direction, but the sequence of wave numbers is simply geometric. 
Replacing k v = A" v and k F = A" F in Eq. (IE. lb gives 

N^ 11 =n v -n F ^lnRe. (E.3) 

Therefore the number of grid points in a shell model is about Re 9/4 less than in a 3D DNS. Note also that the number of 
arithmetic operations per variable is about 10 times smaller for a shell model than for a DNS where finite differences 
are calculated in three directions. 

An estimate of the number of time steps Mud is given by the ratio tf/t v where t F and t v are the turn-over times 
l/ui at the forcing and dissipation scales. Assuming Kolmogorov turbulence ui oc e 1 / 3 / 1 / 3 , we find tf/t v m (l F /l v ) 2 ^ 3 . 
Applying Eq. (IE. Il l we find 

M HD = t F /ty * Re 1/2 , (E.4) 

which is the same for DNS and shell models. 

In MHD, for Pm < 1 the previous estimates hold. However, for Pm > 1 the number of grid points in each direction 
increases by a factor l v /l n ~ Pm 1 / 2 , leading to 

l F /l n ~ Re 3/4 Pm 1/2 , (E.5) 

and then to 

AC D -(W/,,) 3 -Re 9/4 Pm 3 / 2 . (E.6) 
For an isotropic MHD shell model, replacing k v by A"'' and kf by A" F in Eq. ( IE. 5b leads to 

<hd = "v - n F * I In Re + I In Pm. (E.7) 

So for Pm > 1, the number of grid points in a shell model is about Re^Pm'^ 2 less than in a 3D DNS. 

An estimate of the number of time steps Mmhd is now given by the ratio t F / t n where t n is the magnetic dissipation 
time. By definition the latter is given by t„ = (^/rj. With t v = l^/v and again Z„// l; « Pm'^ 2 then t n = t v . From Eq. dE.4t 
we find 

Mmhd = t F /t„ * Re 1/2 , (E.8) 

which again is the same for a DNS and a shell model. Note that the number of time steps is approximately the same 
in HD and MHD, and does not depend on Pm. 

Numerical integration 

A system of ODEs can be integrated with standard numerical methods like Runge-Kutta. However, special care 
must be taken because shell model system of equations is stiff. Indeed as shown before the characteristic times vary 
considerably between large and small scales by a factor Re'^ 2 if Pm < 1, or Re ^Pm 1 ' 3 if Pm > 1. 

One could think of using a constant time-step equal to the smallest characteristic time of the problem, correspond- 
ing here to the dissipation time. However, because the dissipation rate fluctuates strongly, this can lead to a numerical 
effect of negative viscosity. This is why it is much better to use an adaptive time step, and in addition keeping the 
latter at least one order of magnitude smaller than the dissipation time. 

Another possibility is to split the time-step in two parts, each part solving different physics. During the first part 
only the non-linear terms are integrated, leading to a value for U„ (in HD). During the second part the exact solution 
for the dissipative term is calculated, replacing U„ by U„ exp(-vA: 2 Af). This helps to avoid the negative viscosity effect 
pre viously mentione d, but it reduces accuracy to a first order approximation O(At). 



Stepanovl (120021) found another way to increase the accuracy of the method. To explain the method, we write the 



ODEs system for a HD shell model in the following form 

U„(t) = F„(U m (t), U q (t), t) - p n U n (t), n = l,N, (E.9) 
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where p„ = vk\, and F„(U,„(t), U q (t), f) is the quadratic function describing the interactions between shells m and q 
with shell n. Introducing 

V„(f) = !/»(*) rap(p„t). (E.10) 

the system dE.9b becomes 

7„(0 = F„(V m (t) exp(-p„,f), 7,(0 exp(-p,0, exp(/V)- (E.l 1) 

Starting with i/„ at ? = fo we want to calculate U„ at t — to + At. Setting V n (to) = t/„(fo), we calculate V„(fo + Af) using 
Eq. dE.l U . Then using Eq. dE. 101 ) we find t/„(fo + At) = V n (to + At)exp(-p„At). Now when calculating 7„(?o + AO, 
the product of coefficients exp(-p m t) and exp(-p q t) with exp(p„t) in Eq. dE.l U can lead to a significant residual error 
especially at small scales where p/At are much larger than unity. Therefore the numerical integration of Eq. dE.l 11 1 
needs to be elaborated a little. Starting with a 4th order Runge-Kutta method, we calculate the explicit quadratic form 
for F„ in order to avoid exp(p„Af) in the numerical integration of Eq. dE.l 11 1. This leads to the following numerical 
scheme 

K\ = AtF n (U m (to),U q (t ),to), 

K; = AtF n (UJt ) exp(-p,„At/2) + exp(-p,„At/2)K}j2, U q (t ) sxp{-p q At/2) + &xp(-p q At/2)K\/2, t + At/2), 
Kl = AtF„(U m (t ) exp(- Pm At/2) + K;J2, U q (t ) exp(-p,Af/2) + K%/2, t + At/2), (E.12) 
< = AtF n (U m (t ) exp(-p m At) + exp(-p m At/2)K? n , U q (t ) exp(-p q At) + exp(-p q At /2)K], t + At), 



Unite + At) 



U n (t ) exp(-p„A0 + {Kl exp(-p„A0 + 2K 2 n exp(-p n At/2) + 2K 3 „ exp(-p„Af/2) + <)/6 + 0(Af 5 ). (E.l 3) 



The factors exp(-pjAt) and exp(-p,Af/2) with ; € [l,N] can be precalculated only once. The negative viscosity effect 
is eliminated and the accuracy is up to a fifth order approximation 0(At 5 ). 



More complex methods can be used to deal with stiffness. The VODE time-stepping scheme dBrown et al.LI 19981) 



can be used for shell model simulations. The integration method is a variable-coefficient form of Backward Differen- 
tiation Formula methods. It requires computing the Jacobian of ODEs. For a local shell model the Jacobian is a sparse 
matrix. However, for non-local models the Jacobian is a full matrix and the method loses efficiency. 



Statistics 

In turbulence, the accuracy of e.g. spectral laws and scaling exponents is directly related to the statistics. High 
accuracy requires to record data over a sufficiently long time, thus again increasing the cost of simulation. Compared 
to DNS, another advantage of shell models is the possibility to reach accurate statistic s. 



This was well demonstrated in the study of helicity cascade by iLessinnes et al.l (1201 11) . The challenge was to 
calculate the helicity spectrum H(k) resulting from the sum of two quantities H ± {k) = ±Ak~ 2 ^ + Bk~ 5 t 3 , that is zero 
at leading order and non-zero at next order. The ratio \H(k)\/\H ± (k)\ cc Ar 1 was as small as 10~ 5 . A relative error of 
10~ 6 has been taken in the VODE time scheme. In addition the time fluctuations of \H ± (k)\ were larger than those of 
energy by a factor k, requiring again more data. 

In Fig. lE.35l the results are illustrated for different amounts of data £ = QT, where Q is the number of independent 
runs, each run being performed during time T. For a helical forcing (Fig. IE.35l left panel) we clearly see that H(k) 
does converge to a well-defined spectrum with increasing while // ± (A:) spectra do not change. On the right panel of 
Fig. IE.351 H(k) is getting smaller with increasing ( as expected for a non-helical forcing. 

The degree of convergence follows the power law £~'/ 2 in agreement with the following formula 

cr (a) = (E.14) 

derived for a Gaussian process a, where cr^ a ) is the standard deviation due to finite sampling of the mean value of a, 
cr a the standard deviation of a and the number of independent values of a within the sampling. 
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Figure E.35: Spectra of \H + \, \H | and \H\, for different amounts of data measured by the parameter (, and for two different forcings, helical (left) 
and non-helical (right). 
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